#!/usr/bin/python
# a very naive usage of Python language...

# import section
import re, sys, string, math, os, shutil

# this is the test case
from case import *

print "-----------------------------------------------"
print " FORTRAN CUBE-FLU Generator 1.0"
print "-----------------------------------------------"

# FUNCTIONS

# print string "a" to a previously opened output file named "output".
# Printed lines are split in parts having size less than lungh_riga_max=132.
# The strategy of the function is to try to cut when the line is more than
# lungh_riga=62. it can very rarely fail.
def ow(a):
	lungh_riga = 62
	lungh_riga_max = 124
# 	lungh_riga = 4096  # gfortran
#	lungh_riga_max = 4192

	i_dump = 0
	a_dump = ""
	a_dump_old = ""
#       short lines (maybe already cut)
	if len(a) < lungh_riga_max:
		output.write(a)
		return
#       long lines parts from first to last - 1
	for i in range(0,len(a)):
		a_dump = a_dump+a[i]
		i_dump = i_dump + 1
# 		power cannot be splitted, so * split is forbidden...
#		if i_dump > lungh_riga and re.search("\+|\-|\\\\|\*",a[i]):
#		don't work, too
#		if i_dump > lungh_riga and re.search("\+|\-|\|",a[i]):
#		it works!
		if i_dump > lungh_riga_max:
			print "wrong line splitting!! try to solve the problem in ow()!!"
			print a
			exit(1)
                if i_dump > lungh_riga and re.search("\+|\-| |,|\\\\",a[i]):
#			print "check: "+a_dump_old
			if a_dump_old[0:5] == "!$omp": 
				output.write("!$omp "+a_dump+"&\n")
			elif a_dump_old[0:5] == "!$acc": 
				output.write("!$acc "+a_dump+"&\n")
			else:
				output.write(a_dump+"&\n")
			i_dump = 0
			a_dump_old=a_dump
			a_dump = ""
#       long lines: last
	if len(a_dump) > 0:
#		print "check: "+a_dump_old
		if a_dump_old[0:5] == "!$omp": 
	 		output.write("!$omp "+a_dump)
		elif a_dump_old[0:5] == "!$acc": 
	 		output.write("!$acc "+a_dump)
		else:
	 		output.write(a_dump)

# just print line a, even if larger than 132 columns. 
# it will be cut in later file post-processing
def ow_s(a):
	output.write(a)

# calculates the expected memory used by cubeflu (in Gb). assumes
# real(sp) has size 4-bytes and real(dp) 8-bytes
def memory_usage_calc(my_kind,N_dim,N_equ,N_add_var,space_variables):
	n_variables = 1
	for indice in range(1,N_dim+1):
		nt = space_variables[str(indice)]['nt']
		ng = space_variables[str(indice)]['ng']
		n_variables = n_variables * (nt+2*ng)
	if my_kind == 'sp':
		size_real = 4 # not properly exact but...
	if my_kind == 'dp':
		size_real = 8 # not properly exact but...
# N_dim --> grid ; N_equ --> var,var_k,rhs,rhs_k ; N_add_var --> add_var
	memory = size_real*n_variables*(N_dim+4*N_equ+N_add_var)  
	memory = memory / (1024*1024*1024)  #results in Gb
	return memory

# NOTE: never used! To be canceled
# find the kind of myvar: grid_variables, variables, add_variables or sca_variables
# NOTE: never used! To be canceled
def find_var_type(myvar,space_variables,variables,add_variables,sca_variables):
	output = 'undefined'
	for i,var in space_variables.iteritems():
		if i != 'n':  # avoid 'n' entry which defines the number of variables
			if myvar == "g_"+var['symbol']:
				output = 'gri'
	for i,var in variables.iteritems():
		if i != 'n':
			if myvar == var['symbol']:
				output = 'var'
	for i,var in add_variables.iteritems():
		if i != 'n':
			if myvar == var['symbol']:
				output = 'add'
	for i,var in sca_variables.iteritems():
		if i != 'n':
			if myvar == var['symbol']:
				output = 'sca'
	if output == 'undefined':
		print "\nVariable "+myvar+" not defined. Define or correct the symbol\n"
		exit(1)
	return output

# find the index of a certain space_var. this functions allows for
# generical names for spatial variables (other part does not allow)
# e.g. it takes y and returns 2, because y is the second spatial variable
def find_index_space_var(my_symbol, space_variables):
	n_space_var = space_variables['n'] 
	for k in range(1,n_space_var+1):
		symbol = space_variables[str(k)]['symbol']
		if my_symbol == symbol:
			return k
	print "Error! Space variables with symbol "+str(my_symbol)+" does not exist"
	exit(1)

# creates the comma-separated list of indexes, e.g. "i_x,i_y,i_z"
# allows for naming spatial variables as desired
def crea_indici(space_variables):
	n_dim = space_variables['n']
	indici = "i_"
	for indice in range(1,n_dim):
		symbol = space_variables[str(indice)]['symbol']
		indici += symbol+",i_"
	indice = n_dim
	symbol = space_variables[str(indice)]['symbol']
	indici += symbol
	return indici

# creates index list for middle point, e.g. nt_x/2,nt_y/2,nt_z/2
# used only when printing on screen this point value
# allows for naming spatial variables as desired
def crea_indici_half(space_variables):
	n_dim = space_variables['n']
	indici = "nt_"
	for indice in range(1,n_dim):
		symbol = space_variables[str(indice)]['symbol']
		indici += symbol+"/2,nt_"
	indice = n_dim
	symbol = space_variables[str(indice)]['symbol']
	indici += symbol+"/2"
	return indici

# creates index list for middle point, but using explicit array syntax
# to make it more readable for !$acc update directive, 
# e.g. nt_x/2:nt_x/2,nt_y/2:nt_y/2,nt_z/2:nt_z/2
# used only when printing on screen this point value
# allows for naming spatial variables as desired
def crea_indici_half_per_update(space_variables):
	n_dim = space_variables['n']
	indici = "nt_"
	for indice in range(1,n_dim):
		symbol = space_variables[str(indice)]['symbol']
		indici += symbol+"/2:nt_"+symbol+"/2,nt_"
	indice = n_dim
	symbol = space_variables[str(indice)]['symbol']
	indici += symbol+"/2:nt_"+symbol+"/2"
	return indici

# creates indexes needed to match dual-block configuration
# allows for naming spatial variables as desired
def crea_indici_locali(space_variables,block_number,dir,var_or_rhs):
	n_dim = space_variables['n']
	indici = ""
	for indice in range(1,n_dim+1):
		symbol = space_variables[str(indice)]['symbol']
                if indice == dir:
			if var_or_rhs == 'var':
				if block_number == '001':
					indici += "i_"+symbol+"-iend_"+symbol+"+1,"
				if block_number == '002':
					indici += "i_"+symbol+"-istart_"+symbol+"+ng_"+symbol+"+1,"
			if var_or_rhs == 'rhs':
				if block_number == '001':
					indici += "i_"+symbol+"-iend_"+symbol+"+1,"
				if block_number == '002':
					indici += "i_"+symbol+"-istart_"+symbol+"+1,"
		else:
			indici += "i_"+symbol+"-istart_"+symbol+"+1,"
	return indici[0:len(indici)-1]

# creates other indexes needed to match dual-block configuration
# allows for naming spatial variables as desired
def crea_indici_exchange_mpi(space_variables,send_or_recv,mpi_index_var,point_index):
	n_dim = space_variables['n']
	indici = ""
	if send_or_recv == 'send':
		for ind in range(1,n_dim+1):
			symbol = space_variables[str(ind)]['symbol']
			if ind == mpi_index_var:
				indici += 'i_'+symbol+"- point_get+1,"
			else:
				indici += 'i_'+symbol+"-istart_"+symbol+"+1,"
	return indici[0:len(indici)-1]

# creates indexes needed to read/write parallel I/O. since
# only true boundaries must print ghost ng_start and ng_end
# are used, e.g. istart_x-ng_x_start:iend_x+ng_x_end,...
# allows for naming spatial variables as desired
def crea_indici_start_end_ghost_start_end(space_variables):
	n_dim = space_variables['n']
	indici = "istart_"
	for indice in range(1,n_dim):
		symbol = space_variables[str(indice)]['symbol']
		indici += symbol+"-ng_"+symbol+"_start:iend_"+symbol+"+ng_"+symbol+"_end,istart_"
	indice = n_dim
	symbol = space_variables[str(indice)]['symbol']
	indici += symbol+"-ng_"+symbol+"_start:iend_"+symbol+"+ng_"+symbol+"_end"
	return indici

# creates indexes for post-processing. only spatial variables
# not in average_dir are present here. no indexes at all for
# completely periodic runs.
# allows for naming spatial variables as desired
def crea_indici_post(space_variables,average_dir):
	n_dim = space_variables['n']
	indici = ""
	for indice in range(1,n_dim+1):
		symbol = space_variables[str(indice)]['symbol']
		if average_dir[str(indice)] == 'no':
			indici += "i_"+symbol+","
	if len(indici) != 0:
		return indici[0:len(indici)-1]
	else:
		return indici

# counts how many averaged directions occur
def count_avg_dir(space_variables,average_dir):
	n_dim = space_variables['n']
	n_avg = 0
	for indice in range(1,n_dim+1):
        	if average_dir[str(indice)] == 'yes': 
			n_avg += 1
	return n_avg
			
# transform variables adding the correct suffixes:
# g_x ---------> gri_g_x
# ccc ---------> var_ccc_k
# lap ---------> add_lap_k
# gam ---------> sca_gam_k
def transform_variables(stringa,space_variables,variables,add_variables,sca_variables):
        matrix = "("+crea_indici(space_variables)+")"

	n_sp_var = space_variables['n']
	for ind_sp_var in range(1,n_sp_var+1):
		name_sp_var = space_variables[str(ind_sp_var)]['symbol']
		stringa = stringa.replace("g_"+name_sp_var,"gri_g_"+name_sp_var+matrix)

	n_equ = variables['n']
	for ind_var in range(1,n_equ+1):
		name_var = variables[str(ind_var)]['symbol']
		stringa = stringa.replace(name_var,"var_"+name_var+"_k"+matrix)

#       additional and scalar variables have _k suffix, even if not related to Runge-Kutta
	n_add_var = add_variables['n']
	for ind_add_var in range(1,n_add_var+1):
		name_add_var = add_variables[str(ind_add_var)]['symbol']
		stringa = stringa.replace(name_add_var,"add_"+name_add_var+"_k"+matrix)

	for ind_sca_var,sca_var in sca_variables.iteritems():
		name_sca_var = sca_var['symbol']
		stringa = stringa.replace(name_sca_var,"sca_"+name_sca_var+"_k")
	mod_stringa = stringa
	return mod_stringa

# transform variables adding the correct suffixes:
# g_x     ---------> gri_g_x
# ccc     ---------> var_ccc or tot_ccc  (depending on mpi_or_full)
# der     ---------> tot_der             (total var only if mpi)
# ccm     ---------> pos_ccm             (averaged var)
# lap     ---------> add_lap_k
# gam     ---------> sca_gam_k
# compared to transform_variables it considers mpi_or_full and considers pos and tot var
def transform_variables_post(stringa,space_variables,variables,add_variables,sca_variables,pos_variables,tot_variables,average_dir,mpi_or_full):
        matrix = "("+crea_indici(space_variables)+")"
	indici_post = crea_indici_post(space_variables,average_dir)
	if indici_post != "":
	        matrix_pos = "("+indici_post+")"
	else:
	        matrix_pos = ""

	n_sp_var = space_variables['n']
	for ind_sp_var in range(1,n_sp_var+1):
		name_sp_var = space_variables[str(ind_sp_var)]['symbol']
		stringa = stringa.replace("g_"+name_sp_var,"gri_g_"+name_sp_var+matrix)
	n_equ = variables['n']
	for ind_var in range(1,n_equ+1):
		name_var = variables[str(ind_var)]['symbol']
		if mpi_or_full == 'mpi':
			stringa = stringa.replace(name_var,"var_"+name_var+matrix)
		else:
			stringa = stringa.replace(name_var,"tot_"+name_var+matrix)
	if mpi_or_full == 'full':
		n_tot_var = tot_variables['n']
		for ind_tot_var in range(1,n_tot_var+1):
			name_tot_var = tot_variables[str(ind_tot_var)]['symbol']
			stringa = stringa.replace(name_tot_var,"tot_"+name_tot_var+matrix)

#       additional and scalar variables have _k suffix, even if not related to Runge-Kutta
	n_add_var = add_variables['n']
	for ind_add_var in range(1,n_add_var+1):
		name_add_var = add_variables[str(ind_add_var)]['symbol']
		stringa = stringa.replace(name_add_var,"add_"+name_add_var+"_k"+matrix)
	n_pos_var = pos_variables['n']
	for ind_pos_var in range(1,n_pos_var+1):
		name_pos_var = pos_variables[str(ind_pos_var)]['symbol']
		stringa = stringa.replace(name_pos_var,"pos_"+name_pos_var+matrix_pos)
	for ind_sca_var,sca_var in sca_variables.iteritems():
		name_sca_var = sca_var['symbol']
		stringa = stringa.replace(name_sca_var,"sca_"+name_sca_var+"_k")
	mod_stringa = stringa
	return mod_stringa

# transform variables adding the correct suffixes:
# g_x     ---------> gri_g_x
# ccc     ---------> var_ccc or var_ccc_k
# rhs_ccc ---------> rhs_var_ccc or rhs_var_ccc_k
# lap     ---------> add_lap_k
# gam     ---------> sca_gam_k
# compared to transform_variables it considers k_ornot_k and transforms also rhs
def transform_variables_per_bc(stringa,space_variables,variables,add_variables,sca_variables,k_ornot_k,exchange='no',block_numer='',dir=''):
        matrix = "("+crea_indici(space_variables)+")"

	n_sp_var = space_variables['n']
	for ind_sp_var in range(1,n_sp_var+1):
		name_sp_var = space_variables[str(ind_sp_var)]['symbol']
		stringa = stringa.replace("g_"+name_sp_var,"gri_g_"+name_sp_var+matrix)

# con un artificio la prossima fai sia le var che i rhs
	n_equ = variables['n']
	for ind_var in range(1,n_equ+1):
		name_var = variables[str(ind_var)]['symbol']
		stringa = stringa.replace(name_var,"var_"+name_var+k_ornot_k+matrix)
# artificio per compensare le var messe dove era rhs...
		stringa = stringa.replace("rhs_var","rhs")   

#       additional and scalar variables have _k suffix, even if not related to Runge-Kutta
	n_add_var = add_variables['n']
	for ind_add_var in range(1,n_add_var+1):
		name_add_var = add_variables[str(ind_add_var)]['symbol']
		stringa = stringa.replace(name_add_var,"add_"+name_add_var+"_k"+matrix)

	for ind_sca_var,sca_var in sca_variables.iteritems():
		name_sca_var = sca_var['symbol']
		stringa = stringa.replace(name_sca_var,"sca_"+name_sca_var+"_k")
	mod_stringa = stringa

	if exchange == 'exchange_yes':
        	local_matrix_var = crea_indici_locali(space_variables,block_number,dir,'var')+")"
		mod_stringa_temp = re.sub('EXC_VAR_(\d)','exc_var_r_buf_'+str(block_number)+'(\\1,'+local_matrix_var,stringa)
        	local_matrix_rhs = crea_indici_locali(space_variables,block_number,dir,'rhs')+")"
		mod_stringa = re.sub('EXC_RHS_(\d)','exc_rhs_r_buf_'+str(block_number)+'(\\1,'+local_matrix_rhs,mod_stringa_temp)
		stringa = mod_stringa

	return mod_stringa

# print how many enddo are needed to span all the spatial dimensions
def stampa_enddo(output,space_variables):
	n_dim = space_variables['n']
	for i in range(1,n_dim+1):
		ow("enddo\n")

def stampa_do_for_exchange_mpi(output,space_variables,i_dim,point_index,var_or_rhs):
	n_dim = space_variables['n']
	for i in range(n_dim,0,-1):
		symbol = space_variables[str(i)]['symbol']
		if i != i_dim:
			ow("do i_"+symbol+"=istart_"+symbol+",iend_"+symbol+"\n")
		else:
			if var_or_rhs == 'var':
				ow("do i_"+symbol+" = point_get, point_get+ng_"+symbol+"\n")
			if var_or_rhs == 'rhs':
				ow("do i_"+symbol+" = point_get, point_get\n")

def stampa_do_without_i_dim(output,space_variables,i_dim,start_or_end, ghost_andor_boundary):
	n_dim = space_variables['n']
	if ghost_andor_boundary == 'ghost_and_boundary':
		for i in range(n_dim,0,-1):
			symbol = space_variables[str(i)]['symbol']
			if i != i_dim:
#				ow("do i_"+symbol+"=istart_"+symbol+"-ng_"+symbol+",iend_"+symbol+"+ng_"+symbol+"\n")
				ow("do i_"+symbol+"=istart_"+symbol+",iend_"+symbol+"\n")
			else:
				if start_or_end == 'start':
					ow("do i_"+symbol+"=istart_"+symbol+",istart_"+symbol+"-ng_"+symbol+",-1\n")
				elif start_or_end == 'end':
					ow("do i_"+symbol+"=iend_"+symbol+",iend_"+symbol+"+ng_"+symbol+"\n")
				else:
					print "please specify \'start\' or \'end\' as third parameter of stampa_do_without_i_dim"
					exit(1)
	elif ghost_andor_boundary == 'ghost':
		for i in range(n_dim,0,-1):
			symbol = space_variables[str(i)]['symbol']
			if i != i_dim:
#				ow("do i_"+symbol+"=istart_"+symbol+"-ng_"+symbol+",iend_"+symbol+"+ng_"+symbol+"\n")
				ow("do i_"+symbol+"=istart_"+symbol+",iend_"+symbol+"\n")
			else:
				if start_or_end == 'start':
					ow("do i_"+symbol+"=istart_"+symbol+"-1,istart_"+symbol+"-ng_"+symbol+",-1\n")
				elif start_or_end == 'end':
					ow("do i_"+symbol+"=iend_"+symbol+"+1,iend_"+symbol+"+ng_"+symbol+"\n")
				else:
					print "please specify \'start\' or \'end\' as third parameter of stampa_do_without_i_dim"
					exit(1)
#       gestione solo boundary un po' del cavolo pero'...
	elif ghost_andor_boundary == 'boundary':
		for i in range(n_dim,0,-1):
			symbol = space_variables[str(i)]['symbol']
			if i != i_dim:
				ow("do i_"+symbol+"=istart_"+symbol+",iend_"+symbol+"\n")
#				ow("do i_"+symbol+"=istart_"+symbol+"-ng_"+symbol+",iend_"+symbol+"+ng_"+symbol+"\n")
			else:
				if start_or_end == 'start':
					ow("do i_"+symbol+"=istart_"+symbol+",istart_"+symbol+",-1\n")
				elif start_or_end == 'end':
					ow("do i_"+symbol+"=iend_"+symbol+",iend_"+symbol+"\n")
				else:
					print "please specify \'start\' or \'end\' as third parameter of stampa_do_without_i_dim"
					exit(1)
	else:
		print "errore in ghost_andor_boundary"
		exit(1)
		
def select_bc_random(i_dim,n_dim,space_variables):
	if n_dim != i_dim:
		i_symb_1 = n_dim
	else:
		i_symb_1 = n_dim-1
	if i_symb_1-1 != i_dim:
		i_symb_2 = i_symb_1-1
	else:
		i_symb_2 = i_symb_1-2
	if n_dim > 2:
		symb_2 = space_variables[str(i_symb_2)]['symbol']
		symb_1 = space_variables[str(i_symb_1)]['symbol']
	elif n_dim == 2:
		symb_1 = space_variables[str(i_symb_1)]['symbol']
		symb_2 = ''
	elif n_dim == 1:
		symb_1 = ''
		symb_2 = ''
	return (symb_1,symb_2)
	

def stampa_random(output,symb_1,symb_2):
	if symb_1 == '' and symb_2 == '':  # scalar usage
#		ow("ran_ivec = 1\n")  # not suitable for gpu acceleration
		ow("do ran_ivec = 1,1\n")  # suitable for gpu acceleration
	elif symb_2 == '':     # only one dimensional parallelization
		index_1 = "i_"+symb_1+"- istart_"+symb_1+"+ng_"+symb_1+"+1"
		ow("ran_ivec = ("+index_1+")\n")
	else:
		index_1 = "i_"+symb_1+"- istart_"+symb_1+"+ng_"+symb_1#+"+1"
		index_2 = "i_"+symb_2+"- istart_"+symb_2+"+ng_"+symb_2+"+1"
		amp_2 = "iend_"+symb_2+"- istart_"+symb_2+"+2*ng_"+symb_2+"+1"
		ow("ran_ivec = ("+index_2+")+("+index_1+")*("+amp_2+")\n")
	
	ow("if (ran_idum(ran_ivec).le.0) then\n")
	ow("   ran_idum(ran_ivec)=max(-ran_idum(ran_ivec),1)\n")
	ow("   ran_idum2(ran_ivec)=ran_idum(ran_ivec)\n")
	ow("   do ran_j=ran_NTAB+8,1,-1 ! RANDOM cycle start\n")
	ow("      ran_k=ran_idum(ran_ivec)/ran_IQ1\n")
	ow("      ran_idum(ran_ivec)=ran_IA1*(ran_idum(ran_ivec)-ran_k*ran_IQ1)-ran_k*ran_IR1\n")
	ow("      if (ran_idum(ran_ivec).lt.0) ran_idum(ran_ivec)=ran_idum(ran_ivec)+ran_IM1\n")
	ow("      if (ran_j.le.ran_NTAB) ran_iv(ran_j,ran_ivec)=ran_idum(ran_ivec)\n")
	ow("   enddo ! RANDOM cycle end\n")
	ow("   ran_iy(ran_ivec)=ran_iv(1,ran_ivec)\n")
	ow("endif\n")
	ow("ran_k=ran_idum(ran_ivec)/ran_IQ1\n")
	ow("ran_idum(ran_ivec)=ran_IA1*(ran_idum(ran_ivec)-ran_k*ran_IQ1)-ran_k*ran_IR1\n")
	ow("if (ran_idum(ran_ivec).lt.0) ran_idum(ran_ivec)=ran_idum(ran_ivec)+ran_IM1\n")
	ow("ran_k=ran_idum2(ran_ivec)/ran_IQ2\n")
	ow("ran_idum2(ran_ivec)=ran_IA2*(ran_idum2(ran_ivec)-ran_k*ran_IQ2)-ran_k*ran_IR2\n")
	ow("if (ran_idum2(ran_ivec).lt.0) ran_idum2(ran_ivec)=ran_idum2(ran_ivec)+ran_IM2\n")
	ow("ran_j=1+ran_iy(ran_ivec)/ran_NDIV\n")
	ow("ran_iy(ran_ivec)=ran_iv(ran_j,ran_ivec)-ran_idum2(ran_ivec)\n")
	ow("ran_iv(ran_j,ran_ivec)=ran_idum(ran_ivec)\n")
	ow("if(ran_iy(ran_ivec).lt.1)ran_iy(ran_ivec)=ran_iy(ran_ivec)+ran_IMM1\n")
	ow("ran_ran2=min(ran_AM*ran_iy(ran_ivec),ran_RNMX)\n")
	if symb_1 == '' and symb_2 == '':  # scalar usage
		ow("enddo\n")  # suitable for gpu acceleration
		

def stampa_compact_init(output,compact_operators,space_variables,variables,add_variables,sca_variables):
	ow('subroutine compact_init\n')
	ow('use storage \n')
	ow('implicit none \n')
	ow('integer i,j,k,i1,i2,'+crea_indici(space_variables)+' \n')
	ow('integer nabove_appo,nbelow_appo,m_appo,lda_appo \n')
	ow('real(myk) :: rcond_appo \n')
	n_dim = space_variables['n']
	for i_dim in range(1,n_dim+1):
		var = space_variables[str(i_dim)]
		symb_var = var['symbol']
		nt = str(var['nt'])
		ow("real(myk) :: com_appo_"+symb_var+"(Nt_"+symb_var+",Nt_"+symb_var+")\n")
		ow("integer :: ipvt_appo_"+symb_var+"(Nt_"+symb_var+")\n")
		ow("real(myk) :: z_appo_"+symb_var+"(Nt_"+symb_var+")\n")
		ow("real(myk), dimension(:,:), allocatable :: banded_appo_"+symb_var+"(:,:)\n")
# 1. Ciclo sugli operatori compatti
	for indice,oper in compact_operators.iteritems():
		ow("! compact operator "+indice+"\n")
		if oper['generic'] == 'no':
			print "Operatori compatti non generici non implementati"
			exit(1)
        # 2. Ciclo sulle dimensioni spaziali
		for i_dim in range(1,n_dim+1):
			ow("! direction: "+str(i_dim)+"\n")
			var = space_variables[str(i_dim)]
			symb_var = var['symbol']
  	 	        # 3. Ciclo sulle matrici: A e B
			for AorB in ('A','B'):
				ow("! Matrix "+AorB+"\n")
				ow("com_appo_"+symb_var+" = 0.\n")
				num_bc_1 = oper['bc_1']['n']
		       	        # 4.a Ciclo sul bordo superiore della matrice bc_1
				for i_co in range(1,num_bc_1+1):
					coeff_co = oper['bc_1'][str(i_co)]
					n_above = int(coeff_co[AorB]['nabove'])
					n_below = int(coeff_co[AorB]['nbelow'])
					for i in range(-n_below,n_above+1):
						coeff = coeff_co[AorB]['coeff'][str(i)]
						ow("com_appo_"+symb_var+"("+str(i_co)+","+str(i_co)+"+("+str(i)+")) = "+str(coeff)+"\n")
				num_bc_2 = oper['bc_2']['n']
				# 4.b Ciclo sul bordo inferiore della matrice bc_1
				for i_co in range(-num_bc_2+1,1):  # per arrivare a 0
					coeff_co = oper['bc_2'][str(i_co)]
					n_above = int(coeff_co[AorB]['nabove'])
					n_below = int(coeff_co[AorB]['nbelow'])
					for i in range(-n_below,n_above+1):
						coeff = coeff_co[AorB]['coeff'][str(i)]
						indice_x = "Nt_"+symb_var+"+("+str(i_co)+")"
						indice_y = indice_x+"+"+str(i)
						ow("com_appo_"+symb_var+"(Nt_"+symb_var+"+("+str(i_co)+"),Nt_"+symb_var+"+("+str(i_co)+")+("+str(i)+")) = "+str(coeff)+"\n")
				# 4.c Ciclo sulla parte centrale della matrice
				coeff_co = oper['center']
				n_above = int(coeff_co[AorB]['nabove'])
				n_below = int(coeff_co[AorB]['nbelow'])
	        	        ow("do i_"+symb_var+"=1+"+str(num_bc_1)+",Nt_"+symb_var+"-"+str(num_bc_2)+"\n")
				for i in range(-n_below,n_above+1):
					coeff = coeff_co[AorB]['coeff'][str(i)]
					ow("com_appo_"+symb_var+"(i_"+symb_var+",i_"+symb_var+"+("+str(i)+")) = "+str(coeff)+"\n")
	        	        ow("enddo\n")
				# 4.d trasformazione matrice in forma bandata
				if AorB == 'A':
					n_above = int(oper[AorB]['nabovemax'])
					n_below = int(oper[AorB]['nbelowmax'])
					n_above_name = "nabove_"+AorB+"_"+indice+"_"+symb_var
					n_below_name = "nbelow_"+AorB+"_"+indice+"_"+symb_var
					m_name = "m_"+AorB+"_"+indice+"_"+symb_var
					ow("nabove_appo = "+n_above_name+"\n")
					ow("nbelow_appo = "+n_below_name+"\n")
					ow("m_appo = "+ m_name+"\n")
					ow("lda_appo = 2*nbelow_appo+nabove_appo+1\n")
			                ow("allocate(banded_appo_"+symb_var+"(lda_appo,Nt_"+symb_var+"))\n")
					ow("banded_appo_"+symb_var+" = 0.\n")
                	                ow("do j=1,Nt_"+symb_var+"\n")
                	                ow("   i1 = max(1,j-nabove_appo)\n")
                	                ow("   i2 = min(Nt_"+symb_var+",j+nbelow_appo)\n")
#                	                ow("   i1 = max0(1,j-nabove_appo)\n")
#                	                ow("   i2 = min0(Nt_"+symb_var+",j+nbelow_appo)\n")
                	                ow("   do i = i1,i2\n")
                	                ow("      k = i - j + m_appo\n")
                	                ow("      banded_appo_"+symb_var+"(k,j) = com_appo_"+symb_var+"(i,j)\n")
                	                ow("   enddo\n")
                	                ow("enddo\n")
   					ow("call dgbco(banded_appo_"+symb_var+",lda_appo,Nt_"+symb_var+",nbelow_appo,nabove_appo,ipvt_appo_"+symb_var+",rcond_appo,z_appo_"+symb_var+")\n")
					abd_name = AorB+"_banded_"+indice+"_"+symb_var
					ow(abd_name+" = banded_appo_"+symb_var+"\n")
					ow(AorB+"_ipvt_"+indice +"_"+symb_var+" = ipvt_appo_"+symb_var+"\n")
			                ow("deallocate(banded_appo_"+symb_var+")\n")
				if AorB == 'B':
					ow(AorB+"_com_"+indice+"_"+symb_var +"= com_appo_"+symb_var+"\n")
		#		do j = 1, n
		#		   i1 = max0(1, j-mu)
		#		   i2 = min0(n, j+ml)
		#		   do i = i1, i2
		#		      k = i - j + m
		#		      abd(k,j) = a(i,j)
		#		   enddo
		#		enddo
	ow('end subroutine compact_init\n')

def stampa_compact_inversion(output,compact_operators,indice,lhs,rhs,space_variables,variables,add_variables,sca_variables,sp_var):
# extract operator oper variable names
	abd_name = "A_banded_"+indice+"_"+sp_var
	ipvt_name = "A_ipvt_"+indice +"_"+sp_var
	m_name = "m_A_"+indice+"_"+sp_var
	n_below_name = "nbelow_A_"+indice+"_"+sp_var
	oper = compact_operators[indice]
	n_below = int(oper["A"]['nbelowmax'])
	index_sp_var = find_index_space_var(sp_var, space_variables)
	nt = space_variables[str(index_sp_var)]['nt']

	n_below_name_B = "nbelow_B_"+indice+"_"+sp_var
	n_above_name_B = "nabove_B_"+indice+"_"+sp_var

#	com_name = "B_com_"+indice+"_"+sp_var  # senza prefattore

	sp_var_num = find_index_space_var(sp_var, space_variables)
	if space_variables[str(sp_var_num)]['spacing'] == 'uniform':
		prefactor = compact_operators[indice]['prefactor-uniform']
		if prefactor[0:3] == 'der':
			prefactor_final = "(gri_d_"+sp_var+")**"+prefactor[3]
		else:
			prefactor = prefactor+sp_var
			prefactor_final = transform_variables(prefactor,space_variables,variables,add_variables,sca_variables)
#			prefactor_type = find_var_type(prefactor,space_variables,variables,add_variables,sca_variables)
#			if prefactor_type == 'sca':
#				prefactor_final = prefactor_type +"_"+ prefactor+"_k"
#			else:
#				prefactor_final = prefactor_type +"_"+ prefactor + "_k("+crea_indici(space_variables)+")"
		com_name = "("+prefactor_final+")*B_com_"+indice+"_"+sp_var
	else:
		prefactor = compact_operators[indice]['prefactor']
		if prefactor != '':
			prefactor = prefactor+sp_var
			prefactor_final = transform_variables(prefactor,space_variables,variables,add_variables,sca_variables)
	#		prefactor_type = find_var_type(prefactor,space_variables,variables,add_variables,sca_variables)
	#		if prefactor_type == 'sca':
	#			prefactor_final = prefactor_type +"_"+ prefactor+"_k"
	#		else:
	#			prefactor_final = prefactor_type +"_"+ prefactor + "_k("+crea_indici(space_variables)+")"
	
			com_name = "("+prefactor_final+")*B_com_"+indice+"_"+sp_var
		else:
			com_name = "B_com_"+indice+"_"+sp_var

        ow("do kk=1,Nt_"+sp_var+"\n")
        ow("   "+lhs.replace("i_"+sp_var,"kk")+" = 0.\n")
        ow("   do kb=max(1,kk-"+n_below_name_B+"),min(Nt_"+sp_var+",kk+"+n_above_name_B+")\n")
        ow("      "+lhs.replace("i_"+sp_var,"kk")+" = "+lhs.replace("i_"+sp_var,"kk")+"+"+com_name+"(kk,kb)*"+rhs.replace("i_"+sp_var,"kb")+"\n")
        ow("   enddo\n")
        ow("enddo\n")
#	print"n_below : ----> "+str(n_below)+ "\n"
#	print"nt : ----> "+str(nt)+ "\n"
# first inversion step
	if n_below != 0 and nt > 1 :
		ow("do kk = 1, Nt_"+sp_var+"-1\n")
		ow("   lm = min("+n_below_name+",Nt_"+sp_var+"-kk)  ; l = "+ipvt_name+"(kk)  ; tt = "+lhs.replace("i_"+sp_var,"l")+"\n")
		ow("   if(l .ne. kk) then\n")
		ow("      "+lhs.replace("i_"+sp_var,"l")+" = "+lhs.replace("i_"+sp_var,"kk")+" ;  "+lhs.replace("i_"+sp_var,"kk")+"= tt  \n")
		ow("   endif\n")
		ow("   do ind=1,lm   \n")
		ow("      "+lhs.replace("i_"+sp_var,"kk+ind")+" = "+lhs.replace("i_"+sp_var,"kk+ind")+"+ tt*"+abd_name+"("+m_name+"+ind,kk) \n")
		ow("   enddo  \n")
		ow("enddo  \n")
	ow("do kb = 1, Nt_"+sp_var+"\n")
	ow("   kk = Nt_"+sp_var+"+1-kb ; "+lhs.replace("i_"+sp_var,"kk")+" = "+lhs.replace("i_"+sp_var,"kk")+"/"+abd_name+"("+m_name+",kk) ; lm = min(kk,"+m_name+") -1\n")
	ow("   la = "+m_name+"-lm ; lb = kk-lm ; tt = -"+lhs.replace("i_"+sp_var,"kk")+" ;\n")
	ow("   do ind=1,lm\n")
	ow("      "+lhs.replace("i_"+sp_var,"lb-1+ind")+" = "+lhs.replace("i_"+sp_var,"lb-1+ind")+" + tt*"+abd_name+"(la-1+ind,kk)\n")
	ow("   enddo\n")
	ow("enddo\n")
#nm1 = n - 1
#if (ml .eq. 0) go to 30
#if (nm1 .lt. 1) go to 30
#      do k = 1, nm1
#         lm = min0(ml,n-k)  ; l = ipvt(k)  ; t = b(l)
#         if (l .ne. k) then 
#            b(l) = b(k)  ;  b(k) = t
#         endif
#         do ind=1,lm 
#            b(k+ind) = b(k+ind) + t*abd(m+ind,k)
#         enddo
#      enddo
#30 continue
#!  now solve  u*x = y
#   do kb = 1, n
#      k = n + 1 - kb   ; b(k) = b(k)/abd(m,k)   ; lm = min0(k,m) - 1  ;
#      la = m - lm   ; lb = k - lm   ; t = -b(k)   ;
#      do ind=1,lm
#         b(lb-1+ind) = b(lb-1+ind) + t*abd(la-1+ind,k)
#      enddo
#   enddo
		
# opens the grid file (e.g., grid_x.dat) and calls stampa_grid_dir to print grid 
def stampa_grid(name,space_variables):
	n_dim = space_variables['n']
	for indice in range(1,n_dim+1):
		symbol = space_variables[str(indice)]['symbol']
		output = open(name+"_"+symbol+".dat", 'w')
		stampa_grid_dir(output,indice,space_variables)

# print grid and redefines grid length for non-uniform grids
def stampa_grid_dir(output,indice,space_variables):
		indice_space=indice
		ng = space_variables[str(indice)]['ng']
		symbol = space_variables[str(indice)]['symbol']
		L = space_variables[str(indice)]['L']
		fun = space_variables[str(indice)]['fun'].replace(symbol,'((i-1.)*L/(n_points-1.))')
		n_points = space_variables[str(indice)]['iend']-space_variables[str(indice)]['istart']+1
		for indice in range(1-ng,n_points+ng+1):
			i=indice ; grid = eval(fun)
			output.write(str(grid)+"\n")
			if i== 1:
				starting_grid = grid
			if i== n_points:
				final_grid = grid
# redefining grid length for non-uniform grids
		space_variables[str(indice_space)]['L'] = final_grid-starting_grid


def stampa_read_sol_binary(output,space_variables,variables,sca_variables,visual,mpi_or_not):
	n_dim = space_variables['n']
	ow('subroutine read_sol_binary()\n')
	ow('use storage \n')
        if mpi_or_not == 'mpi_yes':
		ow('use mpi_storage \n')
	ow('implicit none \n')
	indici = crea_indici(space_variables)
	ow("integer "+str(indici)+"\n")
	ow('character(128) nomefile \n')
	ow("nomefile = 'output_???????.binXXXXX' \n")
	ow("write(nomefile(8:14),'(i7.7)') it \n")
        if mpi_or_not == 'mpi_yes':
                ow("write(nomefile(19:23),'(i5.5)') n_rank\n")
	ow("nomefile = trim(path_dati)//trim(nomefile)\n")
	ow("open(unit=10,file=nomefile,form='unformatted')\n")

	n_equ = variables['n']
	ow("   read(10) ")
		
#       print grid
	for indice_space_var in range(1,n_dim+1):
		symbol_space_var = space_variables[str(indice_space_var)]['symbol']
# griglia 1d		ow('gri_g_'+symbol_space_var+"(i_"+symbol_space_var+"),")
		ow('gri_g_'+symbol_space_var+",")
        ow('&\n')
#       print variables
	for indice_var in range(1,n_equ+1):
		symbol_var = variables[str(indice_var)]['symbol']
		if indice_var > 1:
                        ow(", &\n")
		ow("     var_"+symbol_var)
	for ind_sca_var,sca_var in sca_variables.iteritems():
		symbol_sca_var = sca_var['symbol']
		if 'print' in sca_var:
			print_or_not = sca_var['print']
			if print_or_not == 'yes':
				ow(",") # I cannot print only scalar variables so , is mandatory...
				ow("sca_"+symbol_sca_var+"_k")
	ow("\n")
	ow('close(10)\n')
	ow('end subroutine read_sol_binary\n')

def stampa_print_sol_binary(output,space_variables,variables,sca_variables,visual,mpi_or_not):
	n_dim = space_variables['n']
	ow('subroutine print_sol_binary()\n')
	ow('use storage \n')
	if mpi_or_not == 'mpi_yes':
		ow('use mpi_storage \n')
	ow('implicit none \n')
	indici = crea_indici(space_variables)
	ow("integer "+str(indici)+"\n")
	ow('character(128) nomefile \n')
	ow("nomefile = 'output_???????.binXXXXX' \n")
	ow("write(nomefile(8:14),'(i7.7)') it \n")
        if mpi_or_not == 'mpi_yes':
                ow("write(nomefile(19:23),'(i5.5)') n_rank\n")
	ow("nomefile = trim(path_dati)//trim(nomefile)\n")
	ow("open(unit=10,file=nomefile,form='unformatted')\n")
	n_equ = variables['n']
	ow("   write(10) ")
		
#       print grid
	for indice_space_var in range(1,n_dim+1):
		symbol_space_var = space_variables[str(indice_space_var)]['symbol']
# griglia 1d		ow('gri_g_'+symbol_space_var+"(i_"+symbol_space_var+"),")
		ow('gri_g_'+symbol_space_var+",")
        ow('&\n')
#       print variables
	for indice_var in range(1,n_equ+1):
		symbol_var = variables[str(indice_var)]['symbol']
		if indice_var > 1:
                        ow(", &\n")
		ow("      var_"+symbol_var)
	for ind_sca_var,sca_var in sca_variables.iteritems():
		symbol_sca_var = sca_var['symbol']
		if 'print' in sca_var:
			print_or_not = sca_var['print']
			if print_or_not == 'yes':
				ow(",") # I cannot print only scalar variables so , is mandatory...
				ow("sca_"+symbol_sca_var+"_k")
	ow("\n")
	ow('close(10)\n')
	ow('end subroutine print_sol_binary\n')

def stampa_writeread_sol_binary_pio(output,space_variables,variables,sca_variables,visual,mpi_or_not,read_or_write):
	n_dim = space_variables['n']
	if n_dim == 1:
		matrix_empty = "(:)"
	if n_dim == 2:
		matrix_empty = "(:,:)"
	if n_dim == 3:
		matrix_empty = "(:,:,:)"
	n_equ = variables['n']
	ow('subroutine '+read_or_write+'_sol_binary_pio()\n')
	ow('use storage \n')
	if mpi_or_not == 'mpi_yes':
		ow('use mpi_storage \n')
	ow('implicit none \n')
	ow('character(128) nomefile \n')
	indici = crea_indici(space_variables)
	ow("integer "+str(indici)+"\n")
	ow("integer(kind=MPI_OFFSET_KIND) disp \n")
	ow("integer :: gsizes(n_dim) , psizes(n_dim) , lsizes(n_dim) , start_indices(n_dim) \n")
	ow("integer :: thefile, subarray , local_array_size, size_of_reals \n")
	ow("integer(MPI_OFFSET_KIND) :: total_array_size \n")
	ow("real(myk),dimension"+matrix_empty+",allocatable :: buffer\n")
	for ind_sp_var in range(1,n_dim+1):
		symbol_space_var = space_variables[str(ind_sp_var)]['symbol']
		ow("integer :: ng_"+symbol_space_var+"_start,ng_"+symbol_space_var+"_end\n")

	if read_or_write == 'write':
		for ind_sp_var in range(1,n_dim+1):
			symbol_space_var = space_variables[str(ind_sp_var)]['symbol']
			ow("ng_"+symbol_space_var+"_start = 0 ; ng_"+symbol_space_var+"_end = 0\n")
	if read_or_write == 'read':  # each process reads ghost nodes
		for ind_sp_var in range(1,n_dim+1):
			symbol_space_var = space_variables[str(ind_sp_var)]['symbol']
			ow("ng_"+symbol_space_var+"_start = ng_"+symbol_space_var+" ; ng_"+symbol_space_var+"_end = ng_"+symbol_space_var+"\n")

	ow("call MPI_Type_size(mpi_prec, size_of_reals, ierr) \n")

	#ow("size_of_reals = kind(1._myk)\n")
	#if read_or_write == 'write':
	#	print "Warning! Assuming size of reals is equal to Fortran kind"

	local_array_string = "local_array_size = "
	total_array_string = "total_array_size = "
	proc_array_string = ""
	for ind_sp_var in range(1,n_dim+1):
		symbol_space_var = space_variables[str(ind_sp_var)]['symbol']
		ow("gsizes("+str(ind_sp_var)+") = nt_"+symbol_space_var+" + 2*ng_"+symbol_space_var+"\n")
		ow("psizes("+str(ind_sp_var)+") = n_block("+str(ind_sp_var)+")\n")
		ow("if(istart_"+symbol_space_var+" == 1) ng_"+symbol_space_var+"_start = ng_"+symbol_space_var+"\n")
		ow("if(iend_"+symbol_space_var+" == nt_"+symbol_space_var+") ng_"+symbol_space_var+"_end = ng_"+symbol_space_var+"\n")
		ow("lsizes("+str(ind_sp_var)+") = iend_"+symbol_space_var+"-istart_"+symbol_space_var+"+1+ng_"+symbol_space_var+"_start + ng_"+symbol_space_var+"_end\n")
		ow("start_indices("+str(ind_sp_var)+") = istart_"+symbol_space_var+"-1+ng_"+symbol_space_var+"-ng_"+symbol_space_var+"_start \n")
		local_array_string += "lsizes("+str(ind_sp_var)+")*"
		total_array_string += "int(gsizes("+str(ind_sp_var)+"),MPI_OFFSET_KIND)*"
		proc_array_string  += "psizes("+str(ind_sp_var)+")*"
	local_array_string = local_array_string[0:len(local_array_string)-1]+"\n"
	total_array_string = total_array_string[0:len(total_array_string)-1]+"\n"
	proc_array_string  = proc_array_string[0:len(proc_array_string)-1]
	ow(local_array_string)
	ow(total_array_string)
	ow("disp = 0\n")
	
	ow("nomefile = 'output_???????.bin' \n")
	ow("write(nomefile(8:14),'(i7.7)') it \n")
	ow("nomefile = trim(path_dati)//trim(nomefile)\n")

	if read_or_write == 'write':
		mode = "MPI_MODE_WRONLY + MPI_MODE_CREATE"
	if read_or_write == 'read':
		mode = "MPI_MODE_RDONLY"
	ow("call MPI_FILE_OPEN(MPI_COMM_WORLD, nomefile,  "+mode+", & \n")
	ow("                   MPI_INFO_NULL, thefile, ierr) \n")

#       print variables
	matrix = "("+crea_indici_start_end_ghost_start_end(space_variables)+")"
	ow("allocate(buffer"+matrix+") \n")
	for indice_var in range(1,n_equ+1):
		symbol_var = variables[str(indice_var)]['symbol']
		ow("call MPI_Type_create_subarray(n_dim,gsizes,lsizes,start_indices,MPI_ORDER_FORTRAN, mpi_prec, subarray, ierr)\n")
		ow("call MPI_Type_commit(subarray, ierr)\n")
		ow("call MPI_File_set_view(thefile, disp, mpi_prec, subarray, 'NATIVE', MPI_INFO_NULL, ierr)\n")
#		ow("call MPI_FILE_"+read_or_write+"_ALL(thefile,  & \n")
#		ow("var_"+symbol_var+matrix+", & \n")
		if read_or_write == 'write':
			ow("buffer = var_"+symbol_var+matrix+"\n")
		ow("call MPI_FILE_"+read_or_write+"_ALL(thefile, buffer, & \n")
		ow("local_array_size, mpi_prec, MPI_STATUS_IGNORE, ierr) \n")
		if read_or_write == 'read':
			ow("var_"+symbol_var+matrix+" = buffer\n")
		ow("disp = disp + total_array_size*int(size_of_reals,MPI_OFFSET_KIND) ! cast to avoid overflow! \n")

	for ind_sca_var,sca_var in sca_variables.iteritems():
		symbol_sca_var = sca_var['symbol']
		if 'print' in sca_var:
			print_or_not = sca_var['print']
			if print_or_not == 'yes':
				ow(",") # I cannot print only scalar variables so , is mandatory...
				ow("sca_"+symbol_sca_var+"_k")
	ow("\n")
	ow("call MPI_FILE_CLOSE(thefile, ierr) \n")
	ow("deallocate(buffer) \n")

	ow('end subroutine '+read_or_write+'_sol_binary_pio\n')

def stampa_writeread_sol_binary_pio_stream(output,space_variables,variables,sca_variables,visual,mpi_or_not,read_or_write):
	n_dim = space_variables['n']
	n_equ = variables['n']
	ow('subroutine '+read_or_write+'_sol_binary_pio_stream()\n')
	ow('use storage \n')
	ow('use mpi_storage \n')
	ow('implicit none \n')
	ow('character(128) nomefile \n')
	indici = crea_indici(space_variables)
	ow("integer :: "+str(indici)+"\n\n")
	ow("nomefile = 'output_???????.bin' \n")
	ow("write(nomefile(8:14),'(i7.7)') it \n")
	ow("nomefile = trim(path_dati)//trim(nomefile)\n")
	if read_or_write == 'read':
		ow("#ifdef FROM_BIG_ENDIAN\n")
		ow("open(unit=10,file=nomefile,form='unformatted',access='stream',convert='big_endian')\n")
		ow("#elif FROM_LITLLE_ENDIAN\n")
		ow("open(unit=10,file=nomefile,form='unformatted',access='stream',convert='little_endian')\n")
		ow("#else \n")
		ow("open(unit=10,file=nomefile,form='unformatted',access='stream')\n")
		ow("#endif\n\n")
	else:
		ow("open(unit=10,file=nomefile,form='unformatted',access='stream')\n")
        ow("   "+read_or_write+"(10) ")
	for indice_var in range(1,n_equ+1):
		symbol_var = variables[str(indice_var)]['symbol']
		if indice_var > 1:
			ow(",")
		ow("tot_"+symbol_var)
	ow("\n\nclose(10)\n\n")
	ow('end subroutine '+read_or_write+'_sol_binary_pio_stream\n')

def stampa_read_sol(output,space_variables,variables,visual,mpi_or_not):
	n_dim = space_variables['n']
	ow('subroutine read_sol()\n')
	ow('use storage \n')
	if mpi_or_not == 'mpi_yes':
		ow('use mpi_storage \n')
	ow('implicit none \n')
	indici = crea_indici(space_variables)
	ow("integer "+str(indici)+"\n")
	ow('character(128) nomefile \n')
	ow("nomefile = 'output_???????.datXXXXX' \n")
	ow("write(nomefile(8:14),'(i7.7)') it \n")
        if mpi_or_not == 'mpi_yes':
                ow("write(nomefile(19:23),'(i5.5)') n_rank\n")
	ow("nomefile = trim(path_dati)//trim(nomefile)\n")
	ow("open(unit=10,file=nomefile)\n")

        if mpi_or_not == 'mpi_yes':
		ow("if(n_rank.eq.0) then\n ")
	if visual == 'tecplot' and n_dim in (1,2,3):
		ow("   read(10,*)\n ")
		ow("   read(10,*)\n ")
        if mpi_or_not == 'mpi_yes':
		ow("endif\n ")

	for indice in range(n_dim,0,-1):
		symbol = space_variables[str(indice)]['symbol']
		ow("do i_"+symbol+"=istart_"+symbol+",iend_"+symbol+"\n")
	n_equ = variables['n']
	ow("   read(10,*) ")
		
#       print grid
	for indice_space_var in range(1,n_dim+1):
		symbol_space_var = space_variables[str(indice_space_var)]['symbol']
# griglia 1d		ow('gri_g_'+symbol_space_var+"(i_"+symbol_space_var+"),")
		ow('gri_g_'+symbol_space_var+"("+indici+"),")
        ow('&\n')
#       print variables
	for indice_var in range(1,n_equ+1):
		symbol_var = variables[str(indice_var)]['symbol']
		if indice_var > 1:
			ow(", &\n")
		ow("     var_"+symbol_var+"("+crea_indici(space_variables)+")")
	ow("\n")
#       print enddo loops
	for indice in range(1,n_dim+1):
		if visual == 'gnuplot' and n_dim == 2 and indice == 2:
			ow("   read(10,*) \n")
		ow("enddo\n")
	ow('close(10)\n')
	ow('end subroutine read_sol\n')

def stampa_print_sol(output,space_variables,variables,visual,mpi_or_not):
	n_dim = space_variables['n']
	ow('subroutine print_sol()\n')
	ow('use storage \n')
	if mpi_or_not == 'mpi_yes':
		ow('use mpi_storage \n')
	ow('implicit none \n')
	indici = crea_indici(space_variables)
	ow("integer "+str(indici)+"\n")
	ow('character(128) nomefile \n')
	ow("nomefile = 'output_???????.datXXXXX' \n")
	ow("write(nomefile(8:14),'(i7.7)') it \n")
	if mpi_or_not == 'mpi_yes':
		ow("write(nomefile(19:23),'(i5.5)') n_rank\n")
	ow("nomefile = trim(path_dati)//trim(nomefile)\n")
	ow("open(unit=10,file=nomefile)\n")

	if mpi_or_not == 'mpi_yes':
		ow('if(n_rank.eq.0) then \n')
	if visual == 'tecplot' and n_dim in (1,2,3):
		ow("   write(10,'(a)') 'variables = ")
		n_equ = variables['n']
		for ind_sp_var in range(1,n_dim+1):
			name_sp_var = space_variables[str(ind_sp_var)]['name']
			name_sp_var = name_sp_var.replace(' ','_')
			ow("&\n "+name_sp_var)
		for ind_var in range(1,n_equ+1):
			name_var = variables[str(ind_var)]['name']
			name_var = name_var.replace(' ','_')
			ow("&\n "+name_var)
		ow("'\n")
		ow("   write(10,*) 'zone '  ")
		index_name = ('','I','J','K')
		for ind_sp_var in range(1,n_dim+1):
			symbol_sp_var = space_variables[str(ind_sp_var)]['symbol']
			ow(",' "+str(index_name[ind_sp_var])+" = ',nt_"+symbol_sp_var)
		ow("\n")
#	for ind_var in range(1,n_equ+1):
#		name_var = variables[str(ind_var)]['name']
#		name_var = name_var.replace(' ','_')
#		ow(" "+name_var)
#	ow("'\n")
#       print do loops
	if mpi_or_not == 'mpi_yes':
		ow('endif \n')

	for indice in range(n_dim,0,-1):
		symbol = space_variables[str(indice)]['symbol']
		ow("do i_"+symbol+"=istart_"+symbol+",iend_"+symbol+"\n")
	n_equ = variables['n']
	ow("   write(10,'(100(F20.12,1X))') ")
		
#       print grid
	for indice_space_var in range(1,n_dim+1):
		symbol_space_var = space_variables[str(indice_space_var)]['symbol']
# griglia 1d		ow('gri_g_'+symbol_space_var+"(i_"+symbol_space_var+"),")
		ow('gri_g_'+symbol_space_var+"("+indici+"),")
        ow('&\n')
#       print variables
	for indice_var in range(1,n_equ+1):
		symbol_var = variables[str(indice_var)]['symbol']
		if indice_var > 1:
			ow(", &\n")
		ow("     var_"+symbol_var+"("+crea_indici(space_variables)+")")
	ow("\n")
#       print enddo loops
	for indice in range(1,n_dim+1):
		if visual == 'gnuplot' and n_dim == 2 and indice == 2:
			ow("   write(10,*) \n")
		ow("enddo\n")
	ow('close(10)\n')
	ow('end subroutine print_sol\n')

def stampa_read_grid(output,space_variables,tot_or_not):
	n_dim = space_variables['n']
	if tot_or_not == 'tot':
		ow('subroutine input_grid_tot\n')
	else:
		ow('subroutine input_grid\n')
	ow('use storage \n')
	ow('implicit none \n')
	ow('real(myk) :: appo_a, appo_b \n')
	indici = crea_indici(space_variables)
	ow("integer "+str(indici)+"\n")
	for indice in range(1,n_dim+1):
		symbol = space_variables[str(indice)]['symbol']
		n_points = space_variables[str(indice)]['iend']-space_variables[str(indice)]['istart']+1
		ow('real(myk), dimension(1-ng_'+symbol+':nt_'+symbol+'+ng_'+symbol+') :: gri_appo_'+symbol+"\n")
	for indice in range(1,n_dim+1):
		ng = space_variables[str(indice)]['ng']
		symbol = space_variables[str(indice)]['symbol']
		L = space_variables[str(indice)]['L']
		fun = space_variables[str(indice)]['fun'].replace(symbol,'((i-1)*L/(n_points-1))')
		n_points = space_variables[str(indice)]['iend']-space_variables[str(indice)]['istart']+1
		#ow('read(10,*)\n')
                first = "istart_"+symbol+"-ng_"+symbol
                last = "iend_"+symbol+"+ng_"+symbol
		isy = 'i_'+symbol
		ow("open(unit=10,file='grid_"+symbol+".dat')\n")
                ow("read(10,*) appo_a \n")
                ow("read(10,*) appo_b \n")
		ow('close(10)'+'\n')
		ow("gri_d_"+symbol+" = 1./(appo_b - appo_a) ! useful only for uniform grids \n")
		ow("open(unit=10,file='grid_"+symbol+".dat')\n")
		ow("do i_"+symbol+"=1-ng_"+symbol+",Nt_"+symbol+"+ng_"+symbol+"\n")
		if tot_or_not == 'tot':
			ow('   read(10,*) gri_appo_'+symbol+'(i_'+symbol+')\n')
		else:
			ow('   if(('+isy+'.ge.'+first+').and.('+isy+'.le.'+last+')) then\n')
			ow('       read(10,*) gri_appo_'+symbol+'(i_'+symbol+')\n')
			ow('   else \n')
			ow('       read(10,*) \n')
			ow('   endif \n')
		ow('enddo'+'\n')
		ow('close(10)'+'\n')
	if tot_or_not != 'tot':
		for indice in range(1,n_dim+1):
			symbol = space_variables[str(indice)]['symbol']
			for indice_bis in range(1,n_dim+1):
				symbol_bis = space_variables[str(indice_bis)]['symbol']
			#	if symbol_bis != symbol:
				ow("do i_"+symbol_bis+"=istart_"+symbol_bis+"-ng_"+symbol_bis+",iend_"+symbol_bis+"+ng_"+symbol_bis+"\n")
			ow("   gri_g_"+symbol+"("+indici+")=gri_appo_"+symbol+"(i_"+symbol+")"+"\n")
			for indice_bis in range(1,n_dim+1):
				ow('enddo'+'\n')
	else:
		for indice in range(1,n_dim+1):
			symbol = space_variables[str(indice)]['symbol']
			ow("do i_"+symbol+"=1-ng_"+symbol+",nt_"+symbol+"+ng_"+symbol+"\n")
			ow("   tot_g_"+symbol+"(i_"+symbol+")=gri_appo_"+symbol+"(i_"+symbol+")"+"\n")
			ow('enddo'+'\n')
	if tot_or_not == 'tot':
		ow('end subroutine input_grid_tot\n')
	else:
		ow('end subroutine input_grid\n')
	
# Finds the variables of a loop starting from f[i_loop+1] which need to be privatized 
# The rule is that scalar variables in LHS need to be privatized
# Does not work when conditional if statemente are at line starting
def find_private_variable(f,i_loop,loop_number,all_private_list):
	i_loop_start = i_loop+1
	i_loop_temp = i_loop_start
	private_list = "("
	while i_loop_temp-i_loop_start < loop_number:
		i_loop_temp += 1  
		line = f[i_loop_temp].strip()
		lhs = (line.partition('='))[0].strip()
		if lhs[0:4] == 'sca_' and line.count('='):
			if lhs[0:9] not in private_list.strip():
				private_list = private_list+lhs[0:9]+","
			if lhs[0:9] not in all_private_list:
				all_private_list.append(lhs[0:9])
	if len(private_list) > 0:
		private_true_list = private_list[0:len(private_list)-1] #= "" #private_list[len   .rsplit(',')[0]
	private_true_list += ')'
	return private_true_list

# Finds the number of lines of the loop starting from f[i_loop+1]
def find_do_closure(f,i_loop):
	do_balance = 1
	i_loop_temp = i_loop+1
	while do_balance > 0 :
		i_loop_temp += 1  
		if f[i_loop_temp].strip()[0:2] == 'do': 
			do_balance += 1
		if f[i_loop_temp].strip()[0:5] == 'enddo': 
			do_balance -= 1
	return i_loop_temp-i_loop

# Finds the number of do of the loop starting from f[i_loop+1]
def find_do_number(f,i_loop):
	i_loop_temp = i_loop+1
	found = 0
	while found == 0:
		i_loop_temp += 1  
		if f[i_loop_temp].strip()[0:2] != 'do':
			found = 1
	return i_loop_temp-i_loop-1

# Finds loop to be executed sequentially. Here the only rule is
# that such loops contains exactly once ng_* variables. It is assumed
# that only one loop has to be executed serially.
def find_do_sequential(f,i_loop):
	i_loop_temp = i_loop + 1
	found = 0
	while found == 0 and f[i_loop_temp].strip()[0:2] == 'do':
	#	print "riga seq: "+f[i_loop_temp]
		if f[i_loop_temp].count("ng_") == 1:
			found = 1
		i_loop_temp += 1  
	return i_loop_temp-i_loop-1
	
# Trova la chiusa con la parentesi graffa di una stringa che l'abbia appena aperta
# Esempio: asdasas{ase}ase{esda}eee}Dssdgss}g{ggg}  ha chiusura in posizione 25
# Resituisce pure se la stringa e' un operatore semplice ovvero se contiene 
# al suo interno composizione di operatori
def find_parentesis_closure(stringa):
	n_graffe = 1
        semplice = 1
	closure_position = 0
	while stringa != "":
		s1 = stringa.partition("{")
		s2 = stringa.partition("}")
# debug		print "stringa processata:  "+stringa
		if len(s1[0]) < len(s2[0]):
			semplice = 0
			n_graffe += 1
			stringa = s1[2]
			closure_position += len(s1[0])+1
		else:
			n_graffe -= 1
			closure_position += len(s2[0])+1
			if n_graffe == 0:
				return closure_position,semplice
			else:
				stringa = s2[2]
	print "controllo parentesi graffe fallito"

#prova_stringa = "as}{dasd}{Asd}asd{asd}}"
#print "prova_stringa:  "+prova_stringa
#print "stringa chiusa a carattere:   "+str(find_parentesis_closure(prova_stringa))
#exit(1)
		
# applica operatore ad una espressione semplice (ovvero che non contiene
# altri operatori al suo interno 
def appl_operator(espr,oper,variables,add_variables,sca_variables,space_variables,operators):
# debug	print "oper:  "+oper
# debug	print "espr:  "+espr
	line = ""
	for stride,peso in operators[oper[0:2]]['coeff'].iteritems():

		add_line = "+("+peso+")*("+espr+")"
		if stride == "+n":
			stride = '+nt_'+oper[2]  # aggiunte per inserire condizioni periodiche
		if stride == '-n':
			stride = '-nt_'+oper[2]  # aggiunte per inserire condizioni periodiche

		if operators[oper[0:2]]['generic'] == 'yes':
			add_line = add_line.replace("i_"+oper[2],"i_"+oper[2]+stride)
		else:
			print "Non generic operators are still not implemented"
			exit(1)
		line += add_line

	sp_var_num = find_index_space_var(oper[2], space_variables)
	if space_variables[str(sp_var_num)]['spacing'] == 'uniform':
		prefactor = operators[oper[0:2]]['prefactor-uniform']
		if prefactor[0:3] == 'der':
			line = "((gri_d_"+oper[2]+")**"+prefactor[3]+")*("+line+")"
		else:
			prefactor = prefactor+oper[2]
			prefactor = transform_variables(prefactor,space_variables,variables,add_variables,sca_variables)
			line = "("+prefactor+")*("+line+")"
	else:
		prefactor = operators[oper[0:2]]['prefactor']
		if prefactor != '':
			prefactor = prefactor+oper[2]
			prefactor = transform_variables(prefactor,space_variables,variables,add_variables,sca_variables)
			line = "("+prefactor+")*("+line+")"
	return line

def rhs_make(operators,variables,add_variables,sca_variables,space_variables,stringa):

	stringa = transform_variables(stringa,space_variables,variables,add_variables,sca_variables)
	line = transform_expression(stringa,operators,variables,add_variables,sca_variables,space_variables)

	return line

def rhs_make_post(operators,variables,add_variables,sca_variables,space_variables,pos_variables,tot_variables,stringa,average_dir,mpi_or_full):

	stringa = transform_variables_post(stringa,space_variables,variables,add_variables,sca_variables,pos_variables,tot_variables,average_dir,mpi_or_full)
	line = transform_expression(stringa,operators,variables,add_variables,sca_variables,space_variables)

	return line

# trasforma l'espressione fornita dall'utente in un'espressione
# fortran. L'espressione deve prima essere passata con replace
# di stringa che aggiunga le (i) o di quante dimensioni si abbia
# bisogno. La funzione e' ricorsiva perche' se l'operatore
# non e' semplice (contiene al suo interno altri operatori) 
# la funzione chiama se stessa finche' non si arriva ad un
# operatore semplice. 
# N.B. da testare meglio
def transform_expression(stringa,operators,variables,add_variables,sca_variables,space_variables):
	line = ""

	while stringa != "":
# debug		print "-----------------------------------------------"
# debug		print "stringa in itinere:  "+stringa
# debug		print "-----------------------------------------------"

		sub_stringa = stringa.partition("#")
		line += sub_stringa[0]
		if len(sub_stringa[0]) == len(stringa):  #controllo che ci sia almeno un operatore altrimenti esco
			break
	
# debug		print "len(sub_stringa) :  "+str(len(sub_stringa))+"\n"
# debug		print "len(stringa) :  "+str(len(stringa))+"\n"
# debug		print "sub_stringa[0] :  "+sub_stringa[0]+"\n"
# debug		print "sub_stringa[2] :  "+sub_stringa[2]
	
		oper = (sub_stringa[2].split("]"))[0][1:]
# debug		espr = (sub_stringa[2].split("}"))[0][6:]
		fine_operatore, semplice = find_parentesis_closure(sub_stringa[2][6:])
		espr = sub_stringa[2][6:5+fine_operatore]
# debug		print "prima espr,oper:  "+espr+" "+oper
	
		if semplice == 0:
# debug			print "NON semplice"
			espr = transform_expression(espr,operators,variables,add_variables,sca_variables,space_variables)
# debug			print "dopo espr,oper:  "+espr+" "+oper
# debug		print "semplice"
		line = line+"("+appl_operator(espr,oper,variables,add_variables,sca_variables,space_variables,operators)+")"
# debug		stringa = sub_stringa[2].partition("}")[2]
		stringa = sub_stringa[2][6+fine_operatore:]
# debug		print "line : "+line
# debug		exit(1)
	return line

def assegna_nodi_seriale(space_variables):
	n_dim = space_variables['n']
	for indice in range(1,n_dim+1):
		space_variables[str(indice)]['istart'] = 1
		space_variables[str(indice)]['iend'] = space_variables[str(indice)]['nt']

# add _myk kind type parameter to real numbers, except for first three lines in storage
# (identified by path_dati,sp_kind,dp_kind) and format statements (identified by write or read)
# so not completely portable at the moment
def put_kind_constants(file_nokind,file_kind):
	input_file = open(file_nokind, 'r')
	output_file = open(file_kind, 'w')
	f = input_file.readlines()
	for line in f:
		#if re.search("[path_dati|dp_kind|sp_kind]",line):
		if re.search("path_dati|sp_kind|dp_kind|write\(|read\(|format",line):
			line_kind = line
		else: 
			line_kind = re.sub('(\d+\.)([^\w])','\\1'+'_myk'+'\\2',line)
			line_kind = re.sub('(\d+\.[eE][\+\-]\d+)','\\1'+'_myk',line_kind)
			line_kind = re.sub('([^\w])(\.\d+)([^\w])','\\1'+'\\2'+'_myk'+'\\3',line_kind)
			line_kind = re.sub('([^\w])(\.\d+[eE][\+\-]\d+)','\\1'+'\\2'+'_myk',line_kind)
			line_kind = re.sub('(\d+\.\d+)([^\w])','\\1'+'_myk'+'\\2',line_kind)
			line_kind = re.sub('(\d+\.\d+[eE][\+\-]\d+)','\\1'+'_myk',line_kind)
		output_file.write(line_kind)

def put_block_constants(file_noblock,file_block,block):
	input_file = open(file_noblock, 'r')
	output_file = open(file_block, 'w')
	f = input_file.readlines()
	for line in f:
		if re.search("module storage",line):
			line = re.sub("module storage","module storage_"+block,line)
		if re.search("use storage",line):
			line = re.sub("use storage","use storage_"+block,line)
		if re.search("module mpi_storage",line):
			line = re.sub("module mpi_storage","module mpi_storage_"+block,line)
			if not re.search("end module mpi_storage",line):
				line += "\n\nuse storage_exchange\n"
		if re.search("use mpi_storage",line):
			line = re.sub("use mpi_storage","use mpi_storage_"+block,line)
		if line.strip()[0:10] == 'subroutine':
			line = re.sub("subroutine (\w+)","subroutine \\1_"+block,line)
		if line.strip()[0:14] == 'end subroutine':
			line = re.sub("end subroutine (\w+)","end subroutine \\1_"+block,line)
		if re.search("grid_\w.dat",line):
			line = re.sub("grid_(\w).dat",block+"_grid_\\1.dat",line)
		if re.search("sizes.dat",line):
			line = re.sub("sizes.dat",block+"_sizes.dat",line)
		if re.search("scalars.dat",line):
			line = re.sub("scalars.dat",block+"_scalars.dat",line)
		if re.search("namelist",line):
			line = re.sub("scalars","scalars_"+block,line)
		if re.search(r"nml=scalars",line):
			line = re.sub("scalars","scalars_"+block,line)
		line = re.sub("output_\?\?\?\?\?\?\?.bin'","output_???????.bin_"+block+"'",line)
		line = re.sub("output_\?\?\?\?\?\?\?.datXXXXX'","output_???????.datXXXXX_"+block+"'",line)
		line = re.sub("output_\?\?\?\?\?\?\?.binXXXXX'","output_???????.binXXXXX_"+block+"'",line)
		line = re.sub("total_\?\?\?\?\?\?\?_XXXXX.vtr","total_???????_XXXXX_"+block+".vtr",line)
		line = re.sub("total_\?\?\?\?\?\?\?.pvtr","to"+block+"_???????.pvtr",line)
		line = re.sub("_mpi.dat","_mpi.dat_"+block,line)
		line = re.sub("_full.dat","_full.dat_"+block,line)
		output_file.write(line)

def put_block_constants_post_and_post_it(file_noblock,file_block):
	input_file = open(file_noblock, 'r')
	output_file = open(file_block, 'a')
	f = input_file.readlines()
	for line in f:
		if re.search("'grid_\w.dat'",line):
			line = re.sub("'grid_(\w).dat'","trim(adjustl(BLOCK_NUM_STR(2:5)//'grid_\\1.dat'))",line)
		if re.search(" it ",line):
			line = re.sub(" it "," post_it ",line)
		line = re.sub("'output_\?\?\?\?\?\?\?.bin'","'output_???????.bin'//BLOCK_NUM_STR(1:4)",line)
		line = re.sub("'output_\?\?\?\?\?\?\?.datXXXXX'","'output_???????.datXXXXX'//BLOCK_NUM_STR(1:4)",line)
		line = re.sub("'output_\?\?\?\?\?\?\?.binXXXXX'","'output_???????.binXXXXX'//BLOCK_NUM_STR(1:4)",line)
		output_file.write(line)

def stampa_reduce_post(output, space_variables, reduce_dict, mpi_or_full, first_or_second,average_dir,average_type_sub, reduce_or_enlarge):
	ow("subroutine "+reduce_or_enlarge+"_"+mpi_or_full+"_"+average_type_sub+"_"+first_or_second+"()\n\n")
	ow("use mpi_storage \n")
	ow("use storage \n\n")
	n_dim = space_variables['n']
	nt_product_yes = ""
	nt_product_no = ""
	isize_product_yes = ""
	for indice in range(n_dim,0,-1):
		symbol = space_variables[str(indice)]['symbol']
        	if average_dir[str(indice)] == 'yes': 
			nt_product_yes += "nt_"+symbol+"*"
			isize_product_yes += "(iend_"+symbol+"-istart_"+symbol+"+1)*"
		else:
			nt_product_no += "nt_"+symbol+"*"
	if len(nt_product_no) != 0:
		nt_product_no = nt_product_no[0:len(nt_product_no)-1]
	else:
		nt_product_no = str(1)
	n_red = reduce_dict['n']
	for indice in range(1,n_red+1):
		symbol = reduce_dict[str(indice)]

		for index,pos_var in pos_variables.items():
			if type(pos_var) == dict:
				if pos_var['symbol'] == symbol:
					pos_index = index

		average_type = pos_variables[str(pos_index)]['average_type']

		if average_type_sub == average_type:
			if mpi_or_full == 'mpi' and reduce_or_enlarge == 'reduce':
				ow("call MPI_ALLREDUCE(MPI_IN_PLACE,pos_"+symbol+","+nt_product_no+",mpi_prec,MPI_SUM,MPI_COMM_WORLD, ierr)\n")
			if average_type == 'space_and_time':
				if reduce_or_enlarge == 'reduce':
					ow("pos_"+symbol+" = pos_"+symbol+" / ("+nt_product_yes+" ((post_it-post_itmin)/post_itout+1))\n")
				if reduce_or_enlarge == 'enlarge':
					ow("pos_"+symbol+" = pos_"+symbol+" * ("+isize_product_yes+" ((post_it-post_itmin)/post_itout))\n")
			elif average_type == 'space':
				ow("pos_"+symbol+" = pos_"+symbol+" / ("+nt_product_yes+" 1 )\n")
		
	ow("end subroutine "+reduce_or_enlarge+"_"+mpi_or_full+"_"+average_type_sub+"_"+first_or_second+" \n\n")

def stampa_post(output, space_variables, post_dict, mpi_or_full, first_or_second,average_dir,variables,add_variables,sca_variables):
	ow("subroutine post_"+mpi_or_full+"_"+first_or_second+"()\n\n")
	ow("use mpi_storage \n")
	ow("use storage \n\n")
	n_dim = space_variables['n']
        matrix_pos = "("
	for indice in range(n_dim,0,-1):
		symbol = space_variables[str(indice)]['symbol']
		if mpi_or_full == 'mpi':
			ow("do i_"+symbol+"=istart_"+symbol+",iend_"+symbol+"\n")
		else:
			ow("do i_"+symbol+"=1,nt_"+symbol+"\n")
	        if average_dir[str(indice)] == 'no': 
			matrix_pos += "i_"+symbol+","
        matrix_pos += ")"
	n_pos = post_dict['n']
	for indice in range(1,n_pos+1):
		symbol = post_dict[str(indice)]['lhs']
		equation = post_dict[str(indice)]['equation']
		symbol_ind = transform_variables_post(symbol,space_variables,variables,add_variables,sca_variables,pos_variables,tot_variables,average_dir,mpi_or_full)
		ow(symbol_ind+" = "+rhs_make_post(operators,variables,add_variables,sca_variables,space_variables,pos_variables,tot_variables,equation,average_dir,mpi_or_full)+"\n")
	stampa_enddo(output,space_variables)
	ow("end subroutine post_"+mpi_or_full+"_"+first_or_second+" \n\n")



def stampa_total_print_tecplot(output, space_variables, variables, average_dir):
	ow("subroutine total_print_tecplot()\n\n")
	ow("use mpi_storage \n")
	ow("use storage \n\n")
	ow("character(20) :: name_total_file \n ")
	ow("print*,'Starting to print TECPLOT global file iteration:',post_it\n")
	ow("name_total_file = 'total_???????.dat'\n")
	ow("write(name_total_file(7:13),'(I7.7)') post_it \n")
	ow("open(unit=10,file=name_total_file) \n ")
	ow("   write(10,'(a)') 'variables = ")
	n_dim = space_variables['n']
	n_equ = variables['n']
	for ind_sp_var in range(1,n_dim+1):
		name_sp_var = space_variables[str(ind_sp_var)]['name']
		name_sp_var = name_sp_var.replace(' ','_')
		ow(" "+name_sp_var+"&\n")
	for ind_var in range(1,n_equ+1):
		name_var = variables[str(ind_var)]['name']
		name_var = name_var.replace(' ','_')
		ow(" "+name_var+"&\n")
	ow("'\n")
	ow("   write(10,*) 'zone '  ")
	index_name = ('','I','J','K')
	for ind_sp_var in range(1,n_dim+1):
		symbol_sp_var = space_variables[str(ind_sp_var)]['symbol']
		ow(",' "+str(index_name[ind_sp_var])+" = ',nt_"+symbol_sp_var)
	ow("\n")
	for indice in range(n_dim,0,-1):
		symbol = space_variables[str(indice)]['symbol']
		ow("do i_"+symbol+"=1,nt_"+symbol+"\n")
	n_equ = variables['n']
	ow("   write(10,'(100(F20.12,1X))') ")
		
#       print grid
	for indice_space_var in range(1,n_dim+1):
		symbol_space_var = space_variables[str(indice_space_var)]['symbol']
		ow('tot_g_'+symbol_space_var+"(i_"+symbol_space_var+"),")
	ow('& \n')
#       print variables
	for indice_var in range(1,n_equ+1):
		symbol_var = variables[str(indice_var)]['symbol']
		if indice_var > 1:
			ow(",&\n")
		ow("      tot_"+symbol_var+"("+crea_indici(space_variables)+")")
	ow("\n")
#       print enddo loops
	for indice in range(1,n_dim+1):
		ow("enddo\n")
	ow("close(10)\n")
	ow("end subroutine total_print_tecplot\n\n")

def stampa_print_average_tecplot(output, space_variables, pos_variables, average_dir, average_type_sub):
	ow("subroutine print_average_tecplot_"+average_type_sub+"(mpi_or_full)\n\n")
	ow("use mpi_storage \n")
	ow("use storage \n\n")
	ow("integer :: mpi_or_full\n\n")
	ow("character(128) :: name_average_file \n ")
	ow("print*,'Starting to print TECPLOT average file iteration:',post_it\n")
	ow("if(mpi_or_full == 1) name_average_file = 'avg_???????_???????_mpi.dat'\n")
	ow("if(mpi_or_full == 2) name_average_file = 'avg_???????_???????_full.dat'\n")
	if average_type_sub == 'space_and_time':
		ow("write(name_average_file(5:11),'(i7.7)') post_itmin \n")
		ow("write(name_average_file(13:19),'(i7.7)') itmax \n")
	if average_type_sub == 'space':
		ow("write(name_average_file(5:11),'(i7.7)') post_it \n")
		ow("write(name_average_file(13:19),'(i7.7)') post_it \n")
	ow("name_average_file = trim(path_dati)//trim(name_average_file) \n")
	ow("open(unit=101,file=name_average_file) \n")
	ow("write(101,'(a)') 'variables = ")
	n_dim = space_variables['n']
	n_pos = pos_variables['n']
	for ind_sp_var in range(1,n_dim+1):
		if average_dir[str(ind_sp_var)] == 'no':
			name_sp_var = space_variables[str(ind_sp_var)]['name']
			name_sp_var = name_sp_var.replace(' ','_')
			ow(" "+name_sp_var)
	for ind_pos in range(1,n_pos+1):
		if pos_variables[str(ind_pos)]['average_type'] == average_type_sub:
			name_pos = pos_variables[str(ind_pos)]['name']
			name_pos = name_pos.replace(' ','_')
			ow("& \n             &  "+name_pos)
	ow("'\n")
	ow("   write(101,*) 'zone '  ")
	index_name = ('','I','J','K')
	for ind_sp_var in range(1,n_dim+1):
		if average_dir[str(ind_sp_var)] == 'no':
			symbol_sp_var = space_variables[str(ind_sp_var)]['symbol']
			ow(",' "+str(index_name[ind_sp_var])+" = ',nt_"+symbol_sp_var)
	ow("\n")
	for indice in range(n_dim,0,-1):
		if average_dir[str(indice)] == 'no':
			symbol = space_variables[str(indice)]['symbol']
			ow("do i_"+symbol+"=1,nt_"+symbol+"\n")
	n_equ = variables['n']

	if n_pos > 0:
		ow("   write(101,'(100(F20.12,1X))') ")
			
	#       print grid
		grid_no_average_exist = "no"
		stringa = ""
		for indice_space_var in range(1,n_dim+1):
			if average_dir[str(indice_space_var)] == 'no':
				grid_no_average_exist = "yes"
				symbol_space_var = space_variables[str(indice_space_var)]['symbol']
				stringa ='tot_g_'+symbol_space_var+"(i_"+symbol_space_var+"),"
		if grid_no_average_exist == "yes":
			stringa = stringa[0:len(stringa)-1]
		ow(stringa)
	#       print variables
		indici_post = crea_indici_post(space_variables,average_dir)
		if indici_post != "":
		        matrix_pos = "("+indici_post+")"
		else:
		        matrix_pos = ""
                stringa = ""
		first_pos = "yes"
		for indice_var in range(1,n_pos+1):
			symbol_var = pos_variables[str(indice_var)]['symbol']
			if pos_variables[str(indice_var)]['average_type'] == average_type_sub:
				if grid_no_average_exist == 'no' and first_pos == "yes":  # no grid print, all homogeneos
					stringa = "   &\n                  pos_"+symbol_var+matrix_pos
				else:
					stringa = "   &\n                , pos_"+symbol_var+matrix_pos
				first_pos = "no"
				ow(stringa)
		ow("\n")
#       print enddo loops
	for indice in range(1,n_dim+1):
		if average_dir[str(indice)] == 'no':
			ow("enddo\n")
	ow("close(101) \n")
	ow("end subroutine print_average_tecplot_"+average_type_sub+"\n\n")

def stampa_read_average_tecplot(output, space_variables, pos_variables, average_dir, average_type_sub):
	ow("subroutine read_average_tecplot_"+average_type_sub+"(mpi_or_full)\n\n")
	ow("use mpi_storage \n")
	ow("use storage \n\n")
	ow("real(myk) :: temp\n\n")
	ow("integer :: mpi_or_full\n\n")
	ow("character(128) :: name_average_file \n ")
	ow("print*,'Starting to read TECPLOT average file iteration:',post_it\n")
	ow("if(mpi_or_full == 1) name_average_file = 'avg_???????_???????_mpi.dat'\n")
	ow("if(mpi_or_full == 2) name_average_file = 'avg_???????_???????_full.dat'\n")
	if average_type_sub == 'space_and_time':
		ow("write(name_average_file(5:11),'(i7.7)')  post_itmin \n")
		ow("write(name_average_file(13:19),'(i7.7)') itmin \n")
	if average_type_sub == 'space':
		ow("write(name_average_file(5:11),'(i7.7)') post_it \n")
		ow("write(name_average_file(13:19),'(i7.7)') post_it \n")
	ow("name_average_file = trim(path_dati)//trim(name_average_file) \n")
	n_dim = space_variables['n']
	n_pos = pos_variables['n']
	ow("open(unit=101,file=name_average_file) \n")
	ow("read(101,*)\n")
	ow("read(101,*)\n")
	for indice in range(n_dim,0,-1):
		if average_dir[str(indice)] == 'no':
			symbol = space_variables[str(indice)]['symbol']
			ow("do i_"+symbol+"=1,nt_"+symbol+"\n")
	n_equ = variables['n']

	if n_pos > 0:
		ow("   read(101,'(100(F20.12,1X))')  ")
			
	#       print grid
		grid_no_average_exist = "no"
		stringa = ""
		for indice_space_var in range(1,n_dim+1):
			if average_dir[str(indice_space_var)] == 'no':
				grid_no_average_exist = "yes"
				symbol_space_var = space_variables[str(indice_space_var)]['symbol']
				stringa +='temp,'
		if grid_no_average_exist == "yes":
			stringa = stringa[0:len(stringa)-1]
		ow(stringa)
	#       print variables
		indici_post = crea_indici_post(space_variables,average_dir)
		if indici_post != "":
		        matrix_pos = "("+indici_post+")"
		else:
		        matrix_pos = ""
                stringa = ""
                first_pos = "yes"
                for indice_var in range(1,n_pos+1):
                        symbol_var = pos_variables[str(indice_var)]['symbol']
                        if pos_variables[str(indice_var)]['average_type'] == average_type_sub:
                                if grid_no_average_exist == 'no' and first_pos == "yes":  # no grid print, all homogeneos
                                        stringa = "   &\n                  pos_"+symbol_var+matrix_pos
                                else:
                                        stringa = "   &\n                , pos_"+symbol_var+matrix_pos
                                first_pos = "no"
                                ow(stringa)
		ow("\n")
#       print enddo loops
	for indice in range(1,n_dim+1):
		if average_dir[str(indice)] == 'no':
			ow("enddo\n")
	ow("close(101) \n")
	ow("end subroutine read_average_tecplot_"+average_type_sub+"\n\n")

def stampa_total_print_plot3d(output, space_variables, variables, average_dir):
	ow("subroutine total_print_plot3d()\n\n")
	ow("use mpi_storage \n")
	ow("use storage \n\n")
	ow("character(20) :: name_total_file \n ")
	ow("print*,'Starting to print PLOT3D global file iteration:',post_it\n")
	ow("name_total_file = 'total_???????.dat'\n")
	ow("write(name_total_file(7:13),'(I7.7)') post_it \n")
	ow("open(unit=11,file=name_total_file,form='unformatted') \n")
#	ow("rewind (11) \n")

	n_sp_var = space_variables['n']
        n_equ = variables['n']

	symbol = space_variables['1']['symbol']
	ow("write(11) nt_"+symbol)
	for ind_sp_var in range(2,n_sp_var+1):
		symbol = space_variables[str(ind_sp_var)]['symbol']
		ow(",nt_"+symbol)
	ow(","+str(n_equ)+" \n")

	stringa = "write(11) "
        for indice in range(1,n_equ+1):
		for ind_sp_var in range(1,n_sp_var+1):
			stringa += "("
                var_symbol = variables[str(indice)]['symbol']
		stringa += "tot_"+var_symbol+"("+crea_indici(space_variables)+")"
		for ind_sp_var in range(1,n_sp_var+1):
			symbol = space_variables[str(ind_sp_var)]['symbol']
			stringa += ",i_"+symbol+"=1,nt_"+symbol+")"
  		stringa += ","
        stringa = stringa[0:len(stringa)-1]
	ow(stringa+"\n\n")
        ow("close  (11) \n\n")

        ow("print *,'Storing grid' \n")
        ow("open(11,file='grid.bin',form='unformatted') \n")
#        ow("rewind (11) \n")
	symbol = space_variables['1']['symbol']
	ow("write(11) nt_"+symbol)
	for ind_sp_var in range(2,n_sp_var+1):
		symbol = space_variables[str(ind_sp_var)]['symbol']
		ow(",nt_"+symbol)
	ow("\n")

	stringa = "write(11) "
        for indice in range(1,n_sp_var+1):
		for ind_sp_var in range(1,n_sp_var+1):
			stringa += "("
                var_symbol = space_variables[str(indice)]['symbol']
		stringa += "tot_g_"+var_symbol+"(i_"+var_symbol+")"
		for ind_sp_var in range(1,n_sp_var+1):
			symbol = space_variables[str(ind_sp_var)]['symbol']
			stringa += ",i_"+symbol+"=1,nt_"+symbol+")"
  		stringa += ","
        stringa = stringa[0:len(stringa)-1]
	ow(stringa+"\n\n")

	ow("close  (11) \n")
	ow("end subroutine total_print_plot3d\n\n")

def stampa_total_print_vtk(output, space_variables, variables, average_dir):
	ow("subroutine total_print_vtk()\n\n")
	ow("use mpi_storage \n")
	ow("use storage \n")
	ow("use lib_vtk_io \n\n")
	ow("character(20) :: name_total_file\n")
	ow("integer :: E_IO\n")
	ow("\n")
	ow("print*,'Starting to print VTK global file iteration:',post_it\n")
	ow("name_total_file = 'total_???????.vtr'\n")
	ow("write(name_total_file(7:13),'(I7.7)') post_it \n")

	symbol = space_variables['1']['symbol']
	prod_sp_var = "nt_"+symbol
	list_sp_var = "nx1=1,nx2=nt_"+symbol
	shap_sp_var = "1:nt_"+symbol
	list_grid = "X=tot_g_"+symbol+"(1:nt_x)"
	n_sp_var = space_variables['n']
	for ind_sp_var in range(2,n_sp_var+1):
		symbol = space_variables[str(ind_sp_var)]['symbol']
		prod_sp_var += "*nt_"+symbol
		if ind_sp_var == 2:
			list_sp_var += ",ny1=1,ny2=nt_"+symbol
			list_grid += ",Y=tot_g_"+symbol+"(1:nt_y)"
		if ind_sp_var == 3:
			list_sp_var += ",nz1=1,nz2=nt_"+symbol
			list_grid += ",Z=tot_g_"+symbol+"(1:nt_z)"
		shap_sp_var += ",1:nt_"+symbol
	if n_sp_var == 1:
		list_sp_var += ",ny1=1,ny2=1,nz1=1,nz2=1"
		list_grid += ",Y=(/0./),Z=(/0./)"
	if n_sp_var == 2:
		list_sp_var += ",nz1=1,nz2=1"
		list_grid += ",Z=(/0./)"

	ow("#ifdef IBM\n")
	ow("E_IO = VTK_INI_XML(output_format='binary', filename=name_total_file, mesh_topology='RectilinearGrid',"+list_sp_var+")\n")
	ow("#else\n")
	ow("E_IO = VTK_INI_XML(output_format='binary', filename=name_total_file, mesh_topology='RectilinearGrid',"+list_sp_var+")\n")
	ow("#endif\n")
	ow("E_IO = VTK_GEO_XML("+list_sp_var+","+list_grid+")\n")
	ow("E_IO = VTK_DAT_XML(var_location='node', var_block_action='OPEN')\n")
        n_equ = variables['n']
	for indice in range(1,n_equ+1):
		symbol = variables[str(indice)]['symbol']
		ow("E_IO = VTK_VAR_XML(NC_NN="+prod_sp_var+",varname='"+symbol+"',var=reshape(tot_"+symbol+"("+shap_sp_var+"),(/"+prod_sp_var+"/))\n")
	ow("E_IO = VTK_DAT_XML(var_location='node', var_block_action='CLOSE')\n")
	ow("E_IO = VTK_GEO_XML()\n")
	ow("E_IO = VTK_END_XML()\n")

	ow("end subroutine total_print_vtk\n\n")

def stampa_print_vtk(output, space_variables, variables, average_dir,my_kind,solver_or_post):
	if solver_or_post == 'solver':
		iter = 'it'
	if solver_or_post == 'post':
		iter = 'post_it'

	ow("subroutine print_vtk()\n\n")
	ow("! To allow for Paraview parallel reading grid must overlap: \n")
	ow("! hence the usage ot istart_?_vtr possibly different from istart_?\n")
	ow("use mpi_storage \n")
	ow("use storage \n")
	ow("use lib_vtk_io \n\n")
	ow("implicit none\n")
	ow("integer :: istart_x_vtr,istart_y_vtr,istart_z_vtr\n")
	ow("integer :: iend_x_vtr,iend_y_vtr,iend_z_vtr\n") 

	ow("character(128) :: name_total_file,name_assembler_file\n")
	ow("integer :: orank,ostart_x,ostart_y,ostart_z,oend_x,oend_y,oend_z,i_cart_x,i_cart_y,i_cart_z\n")
	ow("integer :: E_IO, size_write, i_rank\n\n")
	ow("real(myk),dimension(:),allocatable :: buffer\n")

	if my_kind == 'sp':
		float_bit_size = str(32)
		tp_float = 'Float32'
	if my_kind == 'dp':
		float_bit_size = str(64)
		tp_float = 'Float64'

	ow("print*,'Starting to print PARALLEL VTK global file iteration:',"+iter+"\n")
	ow("name_total_file = 'total_???????_XXXXX.vtr'\n")
	ow("write(name_total_file(7:13),'(I7.7)') "+iter+"\n")
	ow("write(name_total_file(15:19),'(I5.5)') n_rank\n")

	symbol = space_variables['1']['symbol']
	ow("istart_"+symbol+"_vtr = istart_"+symbol+"\n")
	ow("if(istart_"+symbol+" > 1) istart_"+symbol+"_vtr = istart_"+symbol+"-1\n")
	prod_sp_var = "(iend_"+symbol+"-istart_"+symbol+"_vtr+1)"
	list_nt_var = "nt_"+symbol
	shap_sp_var = "istart_"+symbol+"_vtr:iend_"+symbol
	list_sp_var = "nx1=istart_"+symbol+"_vtr,nx2=iend_"+symbol
	list_grid = "X=tot_g_"+symbol+"(istart_"+symbol+"_vtr:iend_"+symbol+")"
	n_sp_var = space_variables['n']
	for ind_sp_var in range(2,n_sp_var+1):
		symbol = space_variables[str(ind_sp_var)]['symbol']
		ow("istart_"+symbol+"_vtr = istart_"+symbol+"\n")
		ow("if(istart_"+symbol+" > 1) istart_"+symbol+"_vtr = istart_"+symbol+"-1\n")
		prod_sp_var += "*(iend_"+symbol+"-istart_"+symbol+"_vtr+1)"
		list_nt_var += ",nt_"+symbol
		shap_sp_var += ",istart_"+symbol+"_vtr:iend_"+symbol
		if ind_sp_var == 2:
			list_sp_var += ",ny1=istart_"+symbol+"_vtr,ny2=iend_"+symbol
			list_grid += ",Y=tot_g_"+symbol+"(istart_"+symbol+"_vtr:iend_"+symbol+")"
		if ind_sp_var == 3:
			list_sp_var += ",nz1=istart_"+symbol+"_vtr,nz2=iend_"+symbol
			list_grid += ",Z=tot_g_"+symbol+"(istart_"+symbol+"_vtr:iend_"+symbol+")"
	ow("size_write = "+prod_sp_var+"\n")
	if n_sp_var == 1:
		list_sp_var += ",ny1=1,ny2=1,nz1=1,nz2=1"
		list_nt_var += ",1,1"
		list_grid += ",Y=(/0./),Z=(/0./)"
	if n_sp_var == 2:
		list_sp_var += ",nz1=1,nz2=1"
		list_nt_var += ",1"
		list_grid += ",Z=(/0./)"

	ow("allocate(buffer(size_write))\n")
	ow("#ifdef IBM\n")
	ow("E_IO = VTK_INI_XML(output_format='binary', filename=name_total_file, mesh_topology='RectilinearGrid',"+list_sp_var+")\n")
	ow("#else\n")
	ow("E_IO = VTK_INI_XML(output_format='binary', filename=name_total_file, mesh_topology='RectilinearGrid',"+list_sp_var+")\n")
	ow("#endif\n")
	ow("E_IO = VTK_GEO_XML("+list_sp_var+","+list_grid+")\n")
	ow("E_IO = VTK_DAT_XML(var_location='node', var_block_action='OPEN')\n")
        n_equ = variables['n']
	for indice in range(1,n_equ+1):
		symbol = variables[str(indice)]['symbol']
		ow("buffer = reshape(var_"+symbol+"("+shap_sp_var+"),(/size_write/))\n")
		ow("E_IO = VTK_VAR_XML(NC_NN=size_write,varname='"+symbol+"',var=buffer)\n")
		#ow("E_IO = VTK_VAR_XML(NC_NN=size_write,varname='"+symbol+"',var=reshape(var_"+symbol+"("+shap_sp_var+"),(/size_write/)))\n")
	ow("E_IO = VTK_DAT_XML(var_location='node', var_block_action='CLOSE')\n")
	ow("E_IO = VTK_GEO_XML()\n")
	ow("E_IO = VTK_END_XML()\n")

	ow("deallocate(buffer)\n")

	ow("if(n_rank == 0) then \n")

	ow("    name_assembler_file = 'total_???????.pvtr'\n")
	ow("    write(name_assembler_file(7:13),'(I7.7)') "+iter+"\n")
	ow("    E_IO = PVTK_INI_XML(filename = name_assembler_file, mesh_topology = 'PRectilinearGrid', tp='Float32', &\n")
	ow("                        nx1=1,nx2=max(1,nt_x),ny1=1,ny2=max(1,nt_y),nz1=1,nz2=max(1,nt_z))\n")
	ow("    name_total_file = 'total_???????_XXXXX.vtr'\n")
	ow("    write(name_total_file(7:13),'(I7.7)') "+iter+"\n")
	ow("    do i_cart_x=0,n_block(1)-1\n")
	ow("    do i_cart_y=0,n_block(2)-1\n")
	ow("    do i_cart_z=0,n_block(3)-1\n")
	ow("       istart_x_vtr = all_st_x(i_cart_x)\n")
	ow("       if(istart_x_vtr > 1) istart_x_vtr = istart_x_vtr-1\n")
	ow("       iend_x_vtr = all_en_x(i_cart_x)\n")
	ow("       istart_y_vtr = all_st_y(i_cart_y)\n")
	ow("       if(istart_y_vtr > 1) istart_y_vtr = istart_y_vtr-1\n")
	ow("       iend_y_vtr = all_en_y(i_cart_y)\n")
	ow("       istart_z_vtr = all_st_z(i_cart_z)\n")
	ow("       if(istart_z_vtr > 1) istart_z_vtr = istart_z_vtr-1\n")
	ow("       iend_z_vtr = all_en_z(i_cart_z)\n")
	ow("       i_rank = cart_to_rank(i_cart_x,i_cart_y,i_cart_z)\n")
	ow("       write(name_total_file(15:19),'(I5.5)') i_rank\n")
	ow("       E_IO = PVTK_GEO_XML(source=name_total_file,   &\n")
	ow("          nx1=istart_x_vtr,nx2=iend_x_vtr,ny1=istart_y_vtr,ny2=iend_y_vtr,nz1=istart_z_vtr,nz2=iend_z_vtr)\n")
	ow("    enddo\n")
	ow("    enddo\n")
	ow("    enddo\n")

	ow("    E_IO = PVTK_DAT_XML(var_location = 'node', var_block_action = 'OPEN')\n")
	for indice in range(1,n_equ+1):
		symbol = variables[str(indice)]['symbol']
		ow("    E_IO = PVTK_VAR_XML(varname = '"+symbol+"', tp='"+tp_float+"')\n")
	ow("    E_IO = PVTK_DAT_XML(var_location = 'node', var_block_action = 'CLOSE')\n")
	ow("    E_IO = PVTK_END_XML()\n")

	ow("endif \n")

	ow("end subroutine print_vtk\n\n")
	
def stampa_storage_allocatable(output,i_kutta,my_kind,space_variables,variables,add_variables,sca_variables,compact_operators,n_decomp,mpi_periodic,all_private_variable):
	n_dim = space_variables['n']
	ow_s("module storage\n\n")
	ow_s("use basic_var \n\n")
	ow_s("use mpi_basic_var \n\n")
	ow_s("character(128), parameter :: path_dati = './'\n\n")

	ow("logical,dimension("+str(n_decomp)+")::periodicbc=(/")
	for i_sp_var in range(1,n_dim):
        	if mpi_periodic[str(i_sp_var)] == 'yes': 
			ow(".true.  ,")
		else:
			ow(".false. ,"),
       	if mpi_periodic[str(n_dim)] == 'yes': 
		ow(".true. ")
	else:
		ow(".false."),
        ow("/)\n")

	if n_decomp == 1:
		ow("real(myk), allocatable, dimension(:,:) :: left_s_buf,right_s_buf\n")
		ow("real(myk), allocatable, dimension(:,:) :: left_r_buf,right_r_buf\n")
	if n_decomp == 2:
		ow("real(myk), allocatable, dimension(:,:,:) :: left_s_buf,right_s_buf,top_s_buf,bottom_s_buf\n")
		ow("real(myk), allocatable, dimension(:,:,:) :: left_r_buf,right_r_buf,top_r_buf,bottom_r_buf\n")
#       next two are not really needed...
		ow("real(myk), allocatable, dimension(:,:,:) :: front_s_buf,back_s_buf\n")
		ow("real(myk), allocatable, dimension(:,:,:) :: front_r_buf,back_r_buf\n")
	if n_decomp == 3:
		ow("real(myk), allocatable, dimension(:,:,:,:) :: left_s_buf,right_s_buf,top_s_buf,bottom_s_buf\n")
		ow("real(myk), allocatable, dimension(:,:,:,:) :: left_r_buf,right_r_buf,top_r_buf,bottom_r_buf\n")
		ow("real(myk), allocatable, dimension(:,:,:,:) :: front_s_buf,back_s_buf\n")
		ow("real(myk), allocatable, dimension(:,:,:,:) :: front_r_buf,back_r_buf\n")

	symbol = space_variables['1']['symbol']
	matrix_var = "(:"
	matrix_rhs = "(:"
	for indice in range(1,n_dim+1):
		symbol = space_variables[str(indice)]['symbol']
                ng = str(space_variables[str(indice)]['ng'])
		nt = str(space_variables[str(indice)]['nt'])
		ow_s("integer :: istart_"+symbol+"\n")
		ow_s("integer :: iend_"+symbol+"\n")
		ow_s("integer :: isize_"+symbol+"\n")
		ow_s("integer :: ng_"+symbol+"\n")
		ow_s("integer :: nt_"+symbol+"\n")
		if indice > 1:
			matrix_var = matrix_var + ",:" 
			matrix_rhs = matrix_rhs + ",:"
	matrix_var += ")"
	matrix_rhs += ")"

        n_avg = count_avg_dir(space_variables,average_dir)
        matrix_pos = "("
        for indice in range(1,n_dim-n_avg+1):
                matrix_pos += ":,"
        matrix_pos = matrix_pos[0:len(matrix_pos)-1]
        matrix_pos += ")"


#       add non existing components otherwise mpi_decompose cannot compile
	for indice in range(n_dim+1,4):
		if indice == 2:
			symbol = 'y'
		if indice == 3:
			symbol = 'z'
		ow_s("integer :: istart_"+symbol+"\n")
		ow_s("integer :: iend_"+symbol+"\n")
		ow_s("integer :: isize_"+symbol+"\n")
		ow_s("integer :: ng_"+symbol+"\n")
		ow_s("integer :: nt_"+symbol+"\n")

	for indice in range(1,n_dim+1):
		symbol = space_variables[str(indice)]['symbol']
		ow("real(myk), dimension(:), allocatable :: tot_g_"+symbol+"\n")

	for indice in range(1,n_dim+1):
		symbol = space_variables[str(indice)]['symbol']
		length  = str(space_variables[str(indice)]['L'])
		ow_s("real(myk) :: l"+symbol+"\n")
#		ow_s("real(myk) :: "+symbol+"(1-ng_"+symbol+":nt_"+symbol+"+ng_"+symbol+")\n")
		ow_s("real(myk), dimension"+matrix_var+", allocatable :: gri_g_"+symbol+"\n")
		ow_s("real(myk) :: gri_d_"+symbol+"\n")

	ow_s("integer :: isize_max_1\n")
#       define isize_max_2 in any case to avoid compilation barks in mpi_decompose.F90
#	if n_dim > 1:
#		ow_s("integer :: isize_max_2\n")
	ow_s("integer :: isize_max_2\n")
        
	if n_dim == 1:
		symbol = space_variables[str(indice)]['symbol']
		nt = str(space_variables[str(indice)]['nt'])

        n_equ = variables['n']
	for indice in range(1,n_equ+1):
		var_symb = variables[str(indice)]['symbol']
		var_symb_k = variables[str(indice)]['symbol']+"_k"
		ow_s("real(myk), dimension"+matrix_var+", allocatable :: var_"+var_symb+" , var_"+var_symb_k+"\n")
	for indice in range(1,n_equ+1):
		var_symb = variables[str(indice)]['symbol']
		var_symb_k = variables[str(indice)]['symbol']+"_k"
		ow_s("real(myk), dimension"+matrix_rhs+", allocatable :: rhs_"+var_symb+" , rhs_"+var_symb_k+"\n")
        n_add = add_variables['n']
	for indice in range(1,n_add+1):
		var_symb_k = add_variables[str(indice)]['symbol']+"_k"
		ow_s("real(myk), dimension"+matrix_var+", allocatable :: add_"+var_symb_k+"\n")
	for ind_sca_var,sca_var in sca_variables.iteritems():
		var_symb_k = sca_var['symbol']+"_k"
		ow_s("real(myk) :: sca_"+var_symb_k+"\n")
	 	ow_s("namelist/scalars/ sca_"+var_symb_k+"\n")

	ow("#ifdef POST_RUNTIME\n")
        n_pos = pos_variables['n']
	n_avg = count_avg_dir(space_variables,average_dir)
#	print "n_avg:  "+str(n_avg)
	for indice in range(1,n_pos+1):
		var_symb = pos_variables[str(indice)]['symbol']
		if n_avg != n_dim:
			ow("real(myk), dimension"+matrix_pos+", allocatable :: pos_"+var_symb+"\n")
		else:
			ow("real(myk) :: pos_"+var_symb+"\n")
	ow("#endif\n")

#       Operatori compatti
	if compact_operators:  # if it not empty
		ow_s("integer kk,lm,la,lb,l,ind,kb\n")
		ow_s("real(myk) :: tt\n")
	for indice,oper in compact_operators.iteritems():
		ow_s("! compact operator "+indice+"\n")
		for i_dim in range(1,n_dim+1):
			ow_s("! direction: "+str(i_dim)+"\n")
			for AorB in ('A','B'):
				var = space_variables[str(i_dim)]
				symb_var = var['symbol']
				n_above = int(oper[AorB]['nabovemax'])
				n_below = int(oper[AorB]['nbelowmax'])
				n_above_name = "nabove_"+AorB+"_"+indice+"_"+symb_var
				n_below_name = "nbelow_"+AorB+"_"+indice+"_"+symb_var
				m_name = "m_"+AorB+"_"+indice+"_"+symb_var
				lda_name = "lda_"+AorB+"_"+indice+"_"+symb_var
				ow_s("integer, parameter :: "+n_above_name +" = "+str(n_above)+"\n")
				ow_s("integer, parameter :: "+n_below_name +" = "+str(n_below)+"\n")
				if AorB == 'A':
					ow_s("integer, parameter :: "+m_name +"="+ str(n_below+n_above+1)+"\n")
					ow_s("integer, parameter :: "+lda_name +"="+ str(2*n_below+n_above+1)+"\n")
					abd_name = AorB+"_banded_"+indice+"_"+symb_var
			                ow_s("real(myk) :: "+abd_name+"("+lda_name+",Nt_"+symb_var+")\n")
					ow_s("integer :: "+AorB+"_ipvt_"+indice +"_"+symb_var+"(Nt_"+symb_var+")\n")
				if AorB == 'B':
					B_name = AorB+"_com_"+indice+"_"+symb_var
			                ow_s("real(myk) :: "+B_name+"(Nt_"+symb_var+",Nt_"+symb_var+")\n")

#       Random storage
#       dynamic allocation --> allocating everything and Fortran code calculates the maximum sizes
	ow_s("\n! random numbers\n")
	ow_s("integer :: ran_nvec\n")

	ow_s("integer, parameter :: ran_IM1=2147483563\n")
	ow_s("integer, parameter :: ran_IM2=2147483399\n")
	ow_s("integer, parameter :: ran_IMM1=ran_IM1-1\n")
	ow_s("integer, parameter :: ran_IA1=40014\n")
	ow_s("integer, parameter :: ran_IA2=40692\n")
	ow_s("integer, parameter :: ran_IQ1=53668\n")
	ow_s("integer, parameter :: ran_IQ2=52774\n")
	ow_s("integer, parameter :: ran_IR1=12211\n")
	ow_s("integer, parameter :: ran_IR2=3791\n")
	ow_s("integer, parameter :: ran_NTAB=32\n")
	ow_s("integer, parameter :: ran_NDIV=1+ran_IMM1/ran_NTAB\n")
	ow_s("real(myk), parameter :: ran_AM=1./ran_IM1\n")
	ow_s("real(myk), parameter :: ran_EPS=1.2e-7  ! cambiare forse per doppia precisione\n")
	ow_s("real(myk), parameter :: ran_RNMX=1.-ran_EPS\n")
	ow_s("integer :: ran_j , ran_k\n")
	ow_s("integer, allocatable :: ran_idum2(:)\n")
	ow_s("integer, allocatable :: ran_iv(:,:)\n")
	ow_s("integer, allocatable :: ran_iy(:)\n")
	ow_s("integer, allocatable :: ran_idum(:)\n")
	ow_s("integer :: ran_ivec\n")
	ow_s("real(myk) :: ran_ran2\n")

	ow_s("\nend module storage\n")
	output.close

# allocate storage subroutine
def stampa_allocate_storage(output,i_kutta,my_kind,space_variables,variables,add_variables,sca_variables,compact_operators,n_decomp, random_ornot, pos_variables, average_dir):
	n_dim = space_variables['n']
	ow("subroutine allocate_storage \n\n")
        ow("use storage\n\n")

	symbol = space_variables['1']['symbol']
	matrix_var = "(istart_"+symbol+"-ng_"+symbol+":iend_"+symbol+"+ng_"+symbol
	matrix_rhs = "(istart_"+symbol+":iend_"+symbol
	matrix_pos = "("
	for indice in range(1,n_dim+1):
		symbol = space_variables[str(indice)]['symbol']
		if indice > 1:
                        matrix_var = matrix_var + ",istart_"+symbol+"-ng_"+symbol+":iend_"+symbol+"+ng_"+symbol 
                        matrix_rhs = matrix_rhs + ",istart_"+symbol+":iend_"+symbol
        	if average_dir[str(indice)] == 'no': 
			matrix_pos += "1:nt_"+symbol+","  
	matrix_var += ")"
	matrix_rhs += ")"
        if len(matrix_pos) > 1:
                matrix_pos = matrix_pos[0:len(matrix_pos)-1]
                matrix_pos += ")"
        else:
                matrix_pos = ""
	for indice in range(1,n_dim+1):
		symbol = space_variables[str(indice)]['symbol']
		ow("allocate(gri_g_"+symbol+matrix_var+")\n")

	for indice in range(1,n_dim+1):
		symbol = space_variables[str(indice)]['symbol']
		ow("allocate(tot_g_"+symbol+"(1-ng_"+symbol+":nt_"+symbol+"+ng_"+symbol+"))\n")

        n_equ = variables['n']
	for indice in range(1,n_equ+1):
		var_symb = variables[str(indice)]['symbol']
		var_symb_k = variables[str(indice)]['symbol']+"_k"
		ow("allocate(var_"+var_symb+matrix_var+" , var_"+var_symb_k+matrix_var+")\n")
	for indice in range(1,n_equ+1):
		var_symb = variables[str(indice)]['symbol']
		var_symb_k = variables[str(indice)]['symbol']+"_k"
		ow("allocate(rhs_"+var_symb+matrix_rhs+" , rhs_"+var_symb_k+matrix_rhs+")\n")
        n_add = add_variables['n']
	for indice in range(1,n_add+1):
		var_symb_k = add_variables[str(indice)]['symbol']+"_k"
		ow("allocate(add_"+var_symb_k+matrix_var+")\n")

	ow("\n")
	if n_decomp == 1:
		ow("allocate(left_s_buf  (ng_x,"+str(n_equ)+"))\n")
		ow("allocate(right_s_buf (ng_x,"+str(n_equ)+"))\n")
		ow("allocate(left_r_buf  (ng_x,"+str(n_equ)+"))\n")
		ow("allocate(right_r_buf (ng_x,"+str(n_equ)+"))\n")
	if n_decomp == 2:
		ow("allocate(left_s_buf  (ng_x,(iend_y-istart_y+1),"+str(n_equ)+"))\n")
		ow("allocate(right_s_buf (ng_x,(iend_y-istart_y+1),"+str(n_equ)+"))\n")
		ow("allocate(bottom_s_buf((iend_x-istart_x+1),ng_y,"+str(n_equ)+"))\n")
		ow("allocate(top_s_buf   ((iend_x-istart_x+1),ng_y,"+str(n_equ)+"))\n")
		ow("allocate(left_r_buf  (ng_x,(iend_y-istart_y+1),"+str(n_equ)+"))\n")
		ow("allocate(right_r_buf (ng_x,(iend_y-istart_y+1),"+str(n_equ)+"))\n")
		ow("allocate(bottom_r_buf((iend_x-istart_x+1),ng_y,"+str(n_equ)+"))\n")
		ow("allocate(top_r_buf   ((iend_x-istart_x+1),ng_y,"+str(n_equ)+"))\n")
	if n_decomp == 3:
		ow("allocate(left_s_buf  (ng_x,(iend_y-istart_y+1),(iend_z-istart_z+1),"+str(n_equ)+"))\n")
		ow("allocate(right_s_buf (ng_x,(iend_y-istart_y+1),(iend_z-istart_z+1),"+str(n_equ)+"))\n")
		ow("allocate(bottom_s_buf((iend_x-istart_x+1),ng_y,(iend_z-istart_z+1),"+str(n_equ)+"))\n")
		ow("allocate(top_s_buf   ((iend_x-istart_x+1),ng_y,(iend_z-istart_z+1),"+str(n_equ)+"))\n")
		ow("allocate(back_s_buf  ((iend_x-istart_x+1),(iend_y-istart_y+1),ng_z,"+str(n_equ)+"))\n")
		ow("allocate(front_s_buf ((iend_x-istart_x+1),(iend_y-istart_y+1),ng_z,"+str(n_equ)+"))\n")
		ow("allocate(left_r_buf  (ng_x,(iend_y-istart_y+1),(iend_z-istart_z+1),"+str(n_equ)+"))\n")
		ow("allocate(right_r_buf (ng_x,(iend_y-istart_y+1),(iend_z-istart_z+1),"+str(n_equ)+"))\n")
		ow("allocate(bottom_r_buf((iend_x-istart_x+1),ng_y,(iend_z-istart_z+1),"+str(n_equ)+"))\n")
		ow("allocate(top_r_buf   ((iend_x-istart_x+1),ng_y,(iend_z-istart_z+1),"+str(n_equ)+"))\n")
		ow("allocate(back_r_buf  ((iend_x-istart_x+1),(iend_y-istart_y+1),ng_z,"+str(n_equ)+"))\n")
		ow("allocate(front_r_buf ((iend_x-istart_x+1),(iend_y-istart_y+1),ng_z,"+str(n_equ)+"))\n")

	ow("#ifdef POST_RUNTIME\n")
        n_pos = pos_variables['n']
	n_avg = count_avg_dir(space_variables,average_dir)
        if n_avg != n_dim:
		for indice in range(1,n_pos+1):
			var_symb = pos_variables[str(indice)]['symbol']
			ow("allocate(pos_"+var_symb+matrix_pos+")\n")
	ow("#endif\n")

#       Random storage
	if random_ornot == "random_serial" or random_ornot == "random_parallel":
	 	ow("RANDOM = .true.\n")
	else:
	 	ow("RANDOM =.false.\n")

	if random_ornot == "random_serial":  # overwriting what imposed in mpi_decompose
		ran_nvec = 1
  
	ow("if(RANDOM) then\n")
	ow("   allocate(ran_idum2(ran_nvec))\n")
	ow("   allocate(ran_iv(ran_NTAB,ran_nvec))\n")
	ow("   allocate(ran_iy(ran_nvec))\n")
	ow("   allocate(ran_idum(ran_nvec))\n")
	ow("   ran_idum2(:) = 123456789\n")
	ow("   ran_iv(:,:) = 0\n")
	ow("   ran_iy(:) = 0\n")
	ow("endif\n")

	ow("\nend subroutine allocate_storage\n")
	output.close

# allocate storage for post-processor
def stampa_allocate_storage_post(output,i_kutta,my_kind,space_variables,variables,add_variables,sca_variables,compact_operators,n_decomp, random_ornot, pos_variables, tot_variables, average_dir):
#pos_variables now refer to FULL_STAT only (MPI_STAT non implemented as concerns them)
	n_dim = space_variables['n']
	ow("subroutine allocate_storage \n\n")
        ow("use storage\n\n")

	symbol = space_variables['1']['symbol']
	matrix_var = "(istart_"+symbol+"-ng_"+symbol+":iend_"+symbol+"+ng_"+symbol
	matrix_rhs = "(istart_"+symbol+":iend_"+symbol
	matrix_pos = "("
	matrix_tot = "(1-ng_"+symbol+":nt_"+symbol+"+ng_"+symbol
	for indice in range(1,n_dim+1):
		symbol = space_variables[str(indice)]['symbol']
		if indice > 1:
                        matrix_var = matrix_var + ",istart_"+symbol+"-ng_"+symbol+":iend_"+symbol+"+ng_"+symbol 
                        matrix_rhs = matrix_rhs + ",istart_"+symbol+":iend_"+symbol
                        matrix_tot = matrix_tot + ",1-ng_"+symbol+":nt_"+symbol+"+ng_"+symbol
        	if average_dir[str(indice)] == 'no': 
			matrix_pos += "1:nt_"+symbol+","  
			#matrix_pos += "1-ng_"+symbol+":nt_"+symbol+"+ng_"+symbol+"," 
	matrix_var += ")"
	matrix_rhs += ")"
	if len(matrix_pos) > 1:
		matrix_pos = matrix_pos[0:len(matrix_pos)-1]
		matrix_pos += ")"
	else:
		matrix_pos = ""
	matrix_tot += ")"

	for indice in range(1,n_dim+1):
		symbol = space_variables[str(indice)]['symbol']
		ow("    allocate(tot_g_"+symbol+"(1-ng_"+symbol+":nt_"+symbol+"+ng_"+symbol+"))\n")

	ow("if(MPI_STAT > 0) then \n")
	for indice in range(1,n_dim+1):
		symbol = space_variables[str(indice)]['symbol']
		ow("  allocate(gri_g_"+symbol+matrix_var+")\n")
        n_equ = variables['n']
	for indice in range(1,n_equ+1):
		var_symb = variables[str(indice)]['symbol']
		ow("  allocate(var_"+var_symb+matrix_var+")\n")
	ow("endif \n")

	ow("if(FULL_STAT > 0) then \n")
	ow("  if(n_rank == 0) then \n")
	for indice in range(1,n_equ+1):
		var_symb = variables[str(indice)]['symbol']
		ow("    allocate(tot_"+var_symb+matrix_tot+")\n")
        n_tot = tot_variables['n']
	for indice in range(1,n_tot+1):
		var_symb = tot_variables[str(indice)]['symbol']
		ow("    allocate(tot_"+var_symb+matrix_tot+")\n")
	ow("  endif \n")
	ow("endif \n")

        n_pos = pos_variables['n']
	n_avg = count_avg_dir(space_variables,average_dir)
        if n_avg != n_dim:
		for indice in range(1,n_pos+1):
			var_symb = pos_variables[str(indice)]['symbol']
			ow("allocate(pos_"+var_symb+matrix_pos+")\n")

#       dead code...
        n_add = add_variables['n']
	for indice in range(1,n_add+1):
		var_symb_k = add_variables[str(indice)]['symbol']+"_k"
		ow("allocate(add_"+var_symb_k+matrix_var+")\n")
	ow("\n")

	ow("\nend subroutine allocate_storage\n")
	output.close

# initialization of pos variables for average post-processing (space-time or space)
def stampa_init_post(output, pos_variables):
	ow("subroutine init_pos \n\n")
        ow("use storage\n\n")

        n_pos = pos_variables['n']
	for indice in range(1,n_pos+1):
		var_symb = pos_variables[str(indice)]['symbol']
		ow("pos_"+var_symb +" = 0.\n")
	ow("\n")

	ow("\nend subroutine init_pos\n")

	ow("\nsubroutine init_pos_space_average \n\n")

        ow("use storage\n\n")

        n_pos = pos_variables['n']
	for indice in range(1,n_pos+1):
		var_symb = pos_variables[str(indice)]['symbol']
		if pos_variables[str(indice)]['average_type'] == 'space':
			ow("pos_"+var_symb +" = 0.\n")
	ow("\n")

	ow("\nend subroutine init_pos_space_average\n\n")

def stampa_enlarge_post_space_and_time_runtime(output, pos_variables):
	ow("\n#ifdef POST_RUNTIME\n")
	ow("subroutine enlarge_post_space_and_time_runtime()\n\n")

        ow("use storage\n\n")

        n_pos = pos_variables['n']
	for indice in range(1,n_pos+1):
		var_symb = pos_variables[str(indice)]['symbol']
		if pos_variables[str(indice)]['average_type'] == 'space_and_time':
			ow("pos_"+var_symb +" = pos_"+var_symb+"*(it-post_itmin-1)/post_itout\n")
	ow("\n")

	ow("end subroutine enlarge_post_space_and_time_runtime\n")
	ow("#endif\n\n")

# initialization of tot variables  for post-processing
def stampa_init_tot(output, tot_variables):
	ow("subroutine init_tot \n\n")
        ow("use storage\n\n")

        n_tot = tot_variables['n']
	for indice in range(1,n_tot+1):
		var_symb = tot_variables[str(indice)]['symbol']
		ow("tot_"+var_symb +" = 0.\n")
	ow("\n")

	ow("\nend subroutine init_tot\n")

# storage module for post-processor
def stampa_storage_allocatable_post(output,i_kutta,my_kind,space_variables,variables,add_variables,sca_variables,compact_operators,n_decomp,mpi_periodic,pos_variables,tot_variables,average_dir):
	n_dim = space_variables['n']
	ow("module storage\n\n")
	ow("use basic_var \n\n")
	ow("use mpi_basic_var \n\n")
	ow("character(128), parameter :: path_dati = './'\n\n")

	ow("logical,dimension("+str(n_decomp)+")::periodicbc=(/")
	for i_sp_var in range(1,n_dim):
        	if mpi_periodic[str(i_sp_var)] == 'yes': 
			ow(".true.  ,")
		else:
			ow(".false. ,"),
       	if mpi_periodic[str(n_dim)] == 'yes': 
		ow(".true. ")
	else:
		ow(".false."),
        ow("/)\n")

	symbol = space_variables['1']['symbol']
	matrix_var = "(:"
	matrix_rhs = "(:"
	for indice in range(1,n_dim+1):
		symbol = space_variables[str(indice)]['symbol']
                ng = str(space_variables[str(indice)]['ng'])
		nt = str(space_variables[str(indice)]['nt'])
		ow("integer :: istart_"+symbol+"\n")
		ow("integer :: iend_"+symbol+"\n")
		ow("integer :: isize_"+symbol+"\n")
		ow("integer :: ng_"+symbol+"\n")
		ow("integer :: nt_"+symbol+"\n")
		if indice > 1:
			matrix_var = matrix_var + ",:" 
			matrix_rhs = matrix_rhs + ",:"
	matrix_var += ")"
	matrix_rhs += ")"

# add non existing components to avoid compilation problems
	for indice in range(n_dim+1,4):
		if indice == 2:
			symbol = 'y'
		if indice == 3:
			symbol = 'z'
		ow_s("integer :: istart_"+symbol+"\n")
		ow_s("integer :: iend_"+symbol+"\n")
		ow_s("integer :: isize_"+symbol+"\n")
		ow_s("integer :: ng_"+symbol+"\n")
		ow_s("integer :: nt_"+symbol+"\n")

	n_avg = count_avg_dir(space_variables,average_dir)
	matrix_pos = "("
	for indice in range(1,n_dim-n_avg+1):
		matrix_pos += ":,"
	matrix_pos = matrix_pos[0:len(matrix_pos)-1]
	matrix_pos += ")"
	for indice in range(1,n_dim+1):
		symbol = space_variables[str(indice)]['symbol']
		length  = str(space_variables[str(indice)]['L'])
		ow("real(myk) :: l"+symbol+"\n")
		ow("real(myk), dimension"+matrix_var+", allocatable :: gri_g_"+symbol+"\n")
		ow("real(myk) :: gri_d_"+symbol+"\n")
	ow("integer :: isize_max_1\n")
#       defines isize_max_2 to avoid compilation problems in mpi_decompose.F90
#	if n_dim > 1:
#		ow_s("integer :: isize_max_2\n")
	ow_s("integer :: isize_max_2\n")
        
	for indice in range(1,n_dim+1):
		symbol = space_variables[str(indice)]['symbol']
		ow("real(myk), dimension(:), allocatable :: tot_g_"+symbol+"\n")
        n_equ = variables['n']
	for indice in range(1,n_equ+1):
		var_symb = variables[str(indice)]['symbol']
		ow("real(myk), dimension"+matrix_var+", allocatable :: var_"+var_symb+"\n")
	for indice in range(1,n_equ+1):
		var_symb = variables[str(indice)]['symbol']
		ow("real(myk), dimension"+matrix_rhs+", allocatable :: tot_"+var_symb+"\n")
        n_tot = tot_variables['n']
	for indice in range(1,n_tot+1):
		var_symb = tot_variables[str(indice)]['symbol']
		ow("real(myk), dimension"+matrix_var+", allocatable :: tot_"+var_symb+"\n")
        n_pos = pos_variables['n']
	n_avg = count_avg_dir(space_variables,average_dir)
#	print "n_avg:  "+str(n_avg)
	for indice in range(1,n_pos+1):
		var_symb = pos_variables[str(indice)]['symbol']
		if n_avg != n_dim:
			ow("real(myk), dimension"+matrix_pos+", allocatable :: pos_"+var_symb+"\n")
		else:
			ow("real(myk) :: pos_"+var_symb+"\n")
        n_add = add_variables['n']
	for indice in range(1,n_add+1):
		var_symb_k = add_variables[str(indice)]['symbol']+"_k"
		ow("real(myk), dimension"+matrix_var+", allocatable :: add_"+var_symb_k+"\n")
	for ind_sca_var,sca_var in sca_variables.iteritems():
		var_symb_k = sca_var['symbol']+"_k"
		ow("real(myk) :: sca_"+var_symb_k+"\n")

	ow("\nreal(myk) :: k_max,dk,k1,k2,k3,k_mod\n")
	ow("real(myk) :: pi = acos(-1.d0)\n")
	ow("real(myk), dimension(:), allocatable :: k_vec\n")
	ow("complex(myk), dimension"+matrix_var+", allocatable :: fou_vex,fou_vey,fou_vez\n")
	ow("real(myk), dimension(:), allocatable :: k1_vec,k2_vec,k3_vec\n")
	if n_avg != n_dim:
		ow("real(myk), dimension"+matrix_pos[0:len(matrix_pos)-1]+",:), allocatable :: spe_vex,spe_vey,spe_vez,spe_vet\n")
	else:
		ow("real(myk), dimension(:), allocatable :: spe_vex,spe_vey,spe_vez,spe_vet\n")

	ow("\nend module storage\n")
	output.close

# modify code to impose openmp and PGI-accelerator directives
def stampa_modify_per_gpu_acc_openmp(output,output_name,i_kutta,space_variables,variables,add_variables,sca_variables,compact_operators,n_decomp,gpuacc_or_openmp,all_private_list,block_number='',mpi_indexes=''):
	n_dim = space_variables['n']
	if mpi_indexes == 'mpi_only_one':
		if n_dim == 1:
			mpi_indexes = "(:,1:1)"
		if n_dim == 2:
			mpi_indexes = "(:,:,1:1)"
		if n_dim == 3:
			mpi_indexes = "(:,:,:,1:1)"

	if gpuacc_or_openmp == 'gpuacc':
		input_name = re.sub('_acc','',output_name)
	if gpuacc_or_openmp == 'openmp':
		input_name = re.sub('_omp','',output_name)

	input_file = open(input_name,'r')
	f=input_file.readlines()
#	for line_i in f:
#	for i in range(0,len(f)):
	i = 0
	while i < len(f):
		line_i = f[i]
		current_line = input_file.readline()
		line = line_i.strip()
		if line_i:
#                       CHANGE NAMES OF SUBROUTINES
			if line[0:10] == 'subroutine' or line[0:14] == 'end subroutine':
				if line[0:14] == 'end subroutine':
					if gpuacc_or_openmp == 'openmp' and not re.search('exchange',output_name):
						ow("!$omp end parallel \n")
				if block_number == '':
					if re.search("\(",line):
						if gpuacc_or_openmp == 'gpuacc':
							line = line.replace('(','_acc(') 
						if gpuacc_or_openmp == 'openmp':
							line = line.replace('(','_omp(') 
						ow(line+"\n") 
					else:
						if gpuacc_or_openmp == 'gpuacc':
							ow(line+"_acc\n")
						if gpuacc_or_openmp == 'openmp':
							ow(line+"_omp\n")
				else:
					if re.search("\(",line):
						if gpuacc_or_openmp == 'gpuacc':
							line = line.replace(str(block_number)+'(','acc_'+str(block_number)+'(') 
						if gpuacc_or_openmp == 'openmp':
							line = line.replace(str(block_number)+'(','omp_'+str(block_number)+'(') 
						ow(line+"\n") 
					else:
						if gpuacc_or_openmp == 'gpuacc':
							line = line.replace(str(block_number),'acc_'+str(block_number)) 
						if gpuacc_or_openmp == 'openmp':
							line = line.replace(str(block_number),'omp_'+str(block_number)) 
						ow(line+"\n") 
#                       PARALLEL_STANDARD
			elif line[0:19]  == '!PARALLEL_STANDARD':
				i_loop = find_do_closure(f,i)
				private_list = find_private_variable(f,i,i_loop,all_private_list)
				if gpuacc_or_openmp == 'gpuacc':
					ow('!$acc region \n')
					if len(private_list) != 1:
						ow('!$acc do independent private '+private_list+'\n')
					else:
						ow('!$acc do independent\n')
				if gpuacc_or_openmp == 'openmp':
					if len(private_list) != 1:
						ow('!$omp do private '+private_list+'\n')
					else:
						ow('!$omp do \n')
				for i_loop_int in range(0,i_loop):
					i = i+1
					ow(f[i])
				if gpuacc_or_openmp == 'gpuacc':
					ow('!$acc end region \n')
				if gpuacc_or_openmp == 'openmp':
					ow('!$omp enddo\n')

#                       PARALLEL_REORDER
			elif line[0:18]  == '!PARALLEL_REORDER':
				i_start_loop = i
				i_loop = find_do_closure(f,i)
				private_list = find_private_variable(f,i,i_loop,all_private_list)
				i_do = find_do_number(f,i)
				i_seq = find_do_sequential(f,i)
				if gpuacc_or_openmp == 'gpuacc':
					ow('!$acc region \n')
					if len(private_list) != 1:
						ow('!$acc do independent private '+private_list+'\n')
					else:
 						ow('!$acc do independent \n')
				if gpuacc_or_openmp == 'openmp':
					if len(private_list) != 1:
						ow('!$omp do private '+private_list+'\n')
					else:
						ow('!$omp do \n')
			#	print "axax: "+str(i_loop)+" ; "+str(i_do)+" ; "+str(i_seq)
				for i_loop_int in range(0,i_seq-1):
					i = i+1
			#		print "1: "+f[i]
					ow(f[i])
				i = i+1
				for i_loop_int in range(i_seq+1,i_do+1):
					i = i+1
			#		print "2: "+f[i]
					ow(f[i])
				if gpuacc_or_openmp == 'gpuacc':
					ow("!$acc do seq \n")
				#print "3: "+f[i_start_loop + i_seq]
				ow(f[i_start_loop + i_seq])
				for i_loop_int in range(i_do+1,i_loop+1):
					i = i+1
			#		print "4: "+f[i]
					ow(f[i])
				if gpuacc_or_openmp == 'gpuacc':
					ow('!$acc end region \n')
				if gpuacc_or_openmp == 'openmp':
					ow('!$omp enddo\n')
#				exit(1)
#                       PARALLEL_ONLYONE
			elif line[0:24]  == '!PARALLEL_ONLYONE_START':
				if gpuacc_or_openmp == 'gpuacc':
					ow('!$acc region \n')
			#TO BE TESTED WITHOUT THIS		ow('!$acc do seq \n')
				if gpuacc_or_openmp == 'openmp':
					ow('!$omp single \n')
# note: end single must be placed after usage of ran_ran2!
			elif line[0:24]  == '!PARALLEL_ONLYONE_END':
				if gpuacc_or_openmp == 'gpuacc':
					ow('!$acc end region \n')
				if gpuacc_or_openmp == 'openmp':  
					ow('!$omp end single \n')
#                       THE LAST (and only) DECLARATION
			elif line[0:7]  == 'integer' and f[i+1][0:7] != 'integer' and not re.search('exchange',output_name):
				ow(f[i])
				if gpuacc_or_openmp == 'openmp':
					ow("!$omp parallel\n")
#                       BEFORE MPI
			elif line[0:17] == 'call mpi_sendrecv' and f[i-3].strip()[0:17] != 'call mpi_sendrecv':
				if gpuacc_or_openmp == 'openmp':
					ow("!$omp end parallel\n")
				if gpuacc_or_openmp == 'gpuacc':
					ow("!$acc update host(left_s_buf"+mpi_indexes+",right_s_buf"+mpi_indexes+")\n")
					if n_decomp >= 2:
						ow("!$acc update host(top_s_buf"+mpi_indexes+",bottom_s_buf"+mpi_indexes+")\n")
					if n_decomp >= 3:
						ow("!$acc update host(front_s_buf"+mpi_indexes+",back_s_buf"+mpi_indexes+")\n")
				ow(line+"\n")
#                       AFTER MPI
                        elif line[0:17] != 'call mpi_sendrecv' and f[i-3].strip()[0:17] == 'call mpi_sendrecv':
				if gpuacc_or_openmp == 'openmp':
					ow("!$omp parallel\n")
				if gpuacc_or_openmp == 'gpuacc':
					ow("!$acc update device(left_r_buf"+mpi_indexes+",right_r_buf"+mpi_indexes+")\n")
					if n_decomp >= 2:
						ow("!$acc update device(top_r_buf"+mpi_indexes+",bottom_r_buf"+mpi_indexes+")\n")
					if n_decomp >= 3:
						ow("!$acc update device(front_r_buf"+mpi_indexes+",back_r_buf"+mpi_indexes+")\n")
				ow(line+"\n")
#                       MPI SEND EXCHANGE
			elif line[0:14] == 'call mpi_isend' and f[i-3].strip()[0:14] != 'call mpi_isend':
				if gpuacc_or_openmp == 'openmp':
					ow("!$omp end parallel\n")
				if gpuacc_or_openmp == 'gpuacc':
					ow("!$acc update host(exc_var_s_buf_"+str(block_number)+")\n")
				ow(line+"\n")
#                       MPI WAIT EXCHANGE
			elif line[0:13] == 'call mpi_wait' and re.search('recv',line):
				ow(line+"\n")
				if gpuacc_or_openmp == 'gpuacc':
					if re.search('rhs exchange',line):
						ow("!$acc update device(exc_var_r_buf_"+str(block_number)+")\n")
					if re.search('var exchange',line):
						ow("!$acc update device(exc_rhs_r_buf_"+str(block_number)+")\n")
					
			elif line[0:13] == 'if((point_get' and f[i+1].strip()[0:4] != 'call':
				ow(f[i])
				if gpuacc_or_openmp == 'openmp':
					ow("!$omp parallel\n")
			elif line[0:18] == 'call mpi_allreduce':
				if gpuacc_or_openmp == 'openmp':
					ow("!$omp end parallel\n")
				ow(f[i])
				if gpuacc_or_openmp == 'openmp':
					ow("!$omp parallel\n")
			else: 
				ow(line+"\n")
		i = i + 1; 

# modify storage module to enable openmp and PGI-accelerator directives
def stampa_modify_storage_per_gpu_acc_openmp(output,output_name,i_kutta,space_variables,variables,add_variables,sca_variables,compact_operators,n_decomp,random_ornot,gpuacc_or_openmp,all_private_list):

	if gpuacc_or_openmp == 'gpuacc':
		input_name = re.sub('_acc','',output_name)
	if gpuacc_or_openmp == 'openmp':
		input_name = re.sub('_omp','',output_name)
	input_file = open(input_name,'r')
	f=input_file.readlines()

	line_prec = ''
	line = ''
	for line_i in f:
		line_prec = line
		line = line_i.strip()
		if gpuacc_or_openmp == 'gpuacc':
			if line[0:10] == 'end module':
				n_sp_var = space_variables['n']
				for ind_sp_var in range(1,n_sp_var+1):
					name_sp_var = space_variables[str(ind_sp_var)]['symbol']
					ow("!$acc mirror(gri_g_"+name_sp_var+")\n")
				
				n_equ = variables['n']
				for ind_var in range(1,n_equ+1):
					name_var = variables[str(ind_var)]['symbol']
					ow("!$acc mirror(var_"+name_var+")\n")
					ow("!$acc mirror(rhs_"+name_var+")\n")
					ow("!$acc mirror(var_"+name_var+"_k)\n")
					ow("!$acc mirror(rhs_"+name_var+"_k)\n")
				
				n_add_var = add_variables['n']
				for ind_add_var in range(1,n_add_var+1):
					name_add_var = add_variables[str(ind_add_var)]['symbol']
					ow("!$acc mirror(add_"+name_add_var+"_k)\n")
				ow("!$acc mirror(left_s_buf,right_s_buf)\n")
				ow("!$acc mirror(left_r_buf,right_r_buf)\n")
				if n_decomp > 2:
					ow("!$acc mirror(top_s_buf,bottom_s_buf)\n")
					ow("!$acc mirror(top_r_buf,bottom_r_buf)\n")
				if n_decomp == 3:
					ow("!$acc mirror(front_s_buf,back_s_buf)\n")
					ow("!$acc mirror(front_r_buf,back_r_buf)\n")
				if random_ornot == "random_yes":
					ow("!$acc mirror(ran_idum2)\n")
					ow("!$acc mirror(ran_iv)\n")
					ow("!$acc mirror(ran_iy)\n")
					ow("!$acc mirror(ran_idum)\n")
				ow(line+"\n")
			else:
				ow(line+"\n")

		if gpuacc_or_openmp == 'openmp':
			written = 'false'
			for i_private in all_private_list:
				if re.search(i_private,line) and re.search('namelist',line):
					ow("#ifndef _OPENMP\n")
					ow(line+"\n")
					ow("#endif\n")
					written = 'true'
			if written == 'false':
				ow(line+"\n")

# initialization of sizes of the problem
def stampa_size_init(output, space_variables):
	ow("subroutine size_init()\n\n")
	ow("use storage\n")
	ow("use mpi_storage\n\n")
	n_dim = space_variables['n']
	ow("logical :: file_exists\n")
        ow("inquire(FILE='sizes.dat', EXIST=file_exists)\n")
        ow("if(file_exists) then\n")
        ow("   write(*,*) 'Reading size from sizes.dat'\n")
        ow("   open(unit=133,FILE='sizes.dat')\n")
	for indice in range(1,n_dim+1):
		symbol = space_variables[str(indice)]['symbol']
		nt = str(space_variables[str(indice)]['nt'])
		ng = str(space_variables[str(indice)]['ng'])
		length  = str(space_variables[str(indice)]['L'])
		ow("   read(133,*) l"+symbol+"\n")
		ow("   read(133,*) nt_"+symbol+"\n")
		ow("   read(133,*) ng_"+symbol+"\n")
        ow("   close(133)\n")
        ow("else\n")
	for indice in range(1,n_dim+1):
		symbol = space_variables[str(indice)]['symbol']
		nt = str(space_variables[str(indice)]['nt'])
		ng = str(space_variables[str(indice)]['ng'])
		length  = str(space_variables[str(indice)]['L'])
		ow("   l"+symbol+"="+length+"\n")
		ow("   nt_"+symbol+"="+nt+"\n")
		ow("   ng_"+symbol+"="+ng+"\n")
        ow("endif\n")
	for indice in range(n_dim+1,4):
		if indice == 2:
			symbol = 'y'
		if indice == 3:
			symbol = 'z'
		ow("l"+symbol+"= 0._myk\n")
		ow("nt_"+symbol+"=1\n")
		ow("ng_"+symbol+"=1\n")
	ow("end subroutine size_init\n\n")

# read iterations to perform and set Runge-Kutta method
def stampa_time_init(output, i_kutta):
	ow("subroutine time_init()\n\n")
	ow("use basic_var\n")
	ow("use mpi_basic_var\n\n")
	ow("open(unit=1,file='iter.dat')\n")
	ow("read(1,*) itmin   ;  read(1,*) itmax   ;   read(1,*) itout ; read(1,*) dt \n")
	ow("#ifdef POST_RUNTIME\n")
	ow("read(1,*) post_itmin   ;  read(1,*) post_itmax ; read(1,*) post_itout\n")
	ow("#endif\n")
	ow("close(1)\n")
	ow("t=0._myk\n")
	ow("it=itmin\n")
	ow("if(n_rank.eq.0) then\n")
	ow("  print*,'--------------------------------------------'\n")
	ow("  print*,'|         CUBE-FLU 0.1 running...          |'\n")
	ow("  print*,'|   starting from iteration: ',itmin,' | '\n")
	ow("  print*,'--------------------------------------------'\n")
	ow("endif\n")
	if i_kutta == 1:
		ow("s=1\n")
		ow("allocate(alpha(s),beta(s))\n")
		ow("alpha=(/0./)\n")
		ow("beta=(/1./)\n")
	if i_kutta == 2:
		ow("s=2\n")
		ow("allocate(alpha(s),beta(s))\n")
		ow("alpha=(/0.,1./)\n")
		ow("beta=(/0.5,0.5/)\n")
	if i_kutta == 3:
		ow("s=2\n")
		ow("allocate(alpha(s),beta(s))\n")
		ow("alpha=(/0.,0.5/)\n")
		ow("beta=(/0.,1./)\n")
	if i_kutta == 4:
		ow("s=3\n")
		ow("allocate(alpha(s),beta(s))\n")
		ow("alpha=(/0.,1./3.,2./3./)\n")
		ow("beta=(/1./4.,0.,3./4./)\n")
	if i_kutta == 5:
		ow("s=4\n")
		ow("allocate(alpha(s),beta(s))\n")
		ow("alpha=(/0.,0.5,0.5,1./)\n")
		ow("beta=(/1./6.,2./6.,2./6.,1./6./)\n")
	ow("if(GPU) then\n")
	ow("!$acc update device(alpha,beta) \n")
	ow("endif\n")
	ow("end subroutine time_init\n\n")

# random number initialization. safe random seeds generated for paralle processes
def stampa_random_init(output,random_ornot):
	ow("subroutine random_init()\n\n")
	ow("use basic_var\n")
	ow("use mpi_basic_var\n\n")
	ow("use storage\n\n")
	ow("implicit none\n")
	ow("integer :: i_rand, il, ir\n\n")
	ow("if(RANDOM) then\n")
	ow("  do i_rand = 1,size(ran_idum)\n")
	ow("     il = i_rand\n")
	if random_ornot == 'random_parallel':
		ow("     ir = n_rank\n")
	else:
		ow("     ir = 1\n")
	ow("     call psdes_safe_s(il,ir)\n")
	ow("     ran_idum(i_rand) = -abs(il)/10\n")
	ow("  enddo\n")
	ow("  print*,'RANDOM VECTOR initialized'\n")
	ow("endif\n")
	ow("end subroutine random_init\n")

# basic var 
def stampa_basic_var(output, my_kind, space_variables,post_or_not,n_decomp):
	n_dim = space_variables['n']
	if post_or_not == 'post_yes':
		ow("module post_basic_var\n\n")
		ow("integer :: MPI_STAT, FULL_STAT, BLOCK_NUM\n")
		ow("character(len=5) :: BLOCK_NUM_STR\n\n")
		ow("end module post_basic_var\n\n")
		ow("module basic_var\n")
		ow("use post_basic_var\n")
	else:
		ow("module basic_var\n")
	ow("\n")
	ow("integer, parameter :: dp_kind = kind(1.d0)\n")
	ow("integer, parameter :: sp_kind = kind(1.)\n")
        ow("integer, parameter :: myk = "+my_kind+"_kind\n\n")
	ow("integer, parameter :: n_dim = "+str(n_dim)+"\n")
	ow("\n")
	ow("integer :: it,itmin,itmax,itout\n")
	if post_or_not == 'post_no':
		ow("#ifdef POST_RUNTIME\n")
	ow("integer :: post_it,post_itmin,post_itmax,post_itout\n")
	if post_or_not == 'post_no':
		ow("#endif\n")
	ow("\n")
	ow("real(myk) :: t,tk,dt\n")
	ow("\n")
	ow("integer :: s\n")
	ow("\n")
	ow("logical :: GPU\n")
	ow("logical :: RANDOM\n")
	ow("\n")
	ow("real(myk), allocatable, dimension(:) :: alpha\n")
	ow("real(myk), allocatable, dimension(:) :: beta\n")
	ow("!$acc mirror(alpha,beta)\n")
	ow("\n")
	if post_or_not == 'post_yes':
		ow("integer :: FFT_STAT, ns  ! number of spectral components\n")
	ow("end module basic_var\n")
	ow("\n")
	ow("module mpi_basic_var\n")
	ow("use mpi\n")
	ow("use basic_var\n")
	ow("\n")
        if my_kind == 'sp':
                ow("integer :: mpi_prec = MPI_REAL\n")
        if my_kind == 'dp':
                ow("integer :: mpi_prec = MPI_REAL8\n")  # passare a DOUBLE_PRECISION forsE
	ow("integer :: n_rank,n_proc,ierr\n")
	ow("\n")
	ow("integer, parameter :: n_decomp = n_dim\n")
	ow("\n")
	ow("integer :: node_dir\n")
	ow("integer :: NGPUXNODE\n")
	ow("integer :: NCOREXNODE\n")
	ow("real    :: ratio_GPU\n")
	ow("integer :: num_devices,my_device\n")
	ow("integer(8) :: device_kind\n")
	ow("integer :: openmp_activate\n")
	ow("integer :: ind_buf_le, ind_buf_ri\n")
#        ow("integer n_left, n_right\n")
	if n_decomp >= 2:
		ow("integer :: ind_buf_to, ind_buf_bo \n")
#        	ow("integer n_bottom, n_top\n")
	if n_decomp == 3:
		ow("integer :: ind_buf_fr, ind_buf_ba\n")
#        	ow("integer n_back, n_front\n")

	ow("\n")
	ow("end module mpi_basic_var\n")

def stampa_init(output,space_variables,variables,add_variables,sca_variables,initial,tipo):
	indent = "   "
	n_dim = space_variables['n']
        matrix = "("+crea_indici(space_variables)+")"
	output.write("subroutine "+tipo+"_init\n")
	ow("use storage\n")
	ow("implicit none\n")
	indici = crea_indici(space_variables)
	ow("integer "+str(indici)+"\n")
        if tipo == 'rhs':  #only one to be parallelized because called during the time iteration
		ow("!PARALLEL_STANDARD\n")
	for indice in range(n_dim,0,-1):
		symbol = space_variables[str(indice)]['symbol']
		if tipo == 'var' or tipo == 'var_k':
			ow("do i_"+symbol+"=istart_"+symbol+"-ng_"+symbol+",iend_"+symbol+"+ng_"+symbol+"\n")
		else:
			ow("do i_"+symbol+"=istart_"+symbol+",iend_"+symbol+"\n")
	n_equ = variables['n']
	for indice in range(1,n_equ+1):
		var_symb = variables[str(indice)]['symbol']
                if tipo == 'var':
			fun = initial[str(indice)]
			fun = transform_variables_per_bc(fun,space_variables,variables,add_variables,sca_variables,'')
#			fun = transform_variables(fun,space_variables,variables,add_variables,sca_variables)
			ow(indent+"var_"+var_symb+matrix+" = "+fun+"\n")
                if tipo == 'var_k':
			#ow(indent+"var_"+var_symb+"_k"+matrix+" = 0.\n")
# meglio non inizializzare a 0 var_k perche' ai bordi poi prima delle bc se divido per 0...
			ow(indent+"var_"+var_symb+"_k"+matrix+" = var_"+var_symb+matrix+"\n")
                if tipo == 'rhs':
			ow(indent+"rhs_"+var_symb+matrix+" = 0.\n")
                if tipo == 'rhs_k':
			ow(indent+"rhs_"+var_symb+"_k"+matrix+" = 0.\n")
	for indice in range(1,n_dim+1):
		ow("enddo\n")
	output.write("end subroutine "+tipo+"_init\n")

def stampa_eval(output,space_variables,variables,add_variables,tipo,rhs_preamble,ini_preamble,sca_variables):
	indent = "   "
	n_dim = space_variables['n']
        matrix = "("+crea_indici(space_variables)+")"
	if tipo in ('ini','var','rhs_k','post'):
		output.write("subroutine "+tipo+"_eval()\n")
	else:
		output.write("subroutine "+tipo+"_eval(ik)\n")
	ow("use storage\n")
#       CHECK THE NEXT....
	ow("use mpi_storage\n")
	ow('implicit none \n')
	indici = crea_indici(space_variables)
	ow("integer :: ik,"+str(indici)+"\n")
#       print rhs_k preamble, i.e. loops to prepare rhs_k evaluation
        if tipo == 'rhs_k' or tipo == 'ini' or tipo == 'post':
		if tipo == 'rhs_k':
			n_preamble = rhs_preamble['n']
		if tipo == 'ini':
			n_preamble = ini_preamble['n']
		if tipo == 'post':
			n_preamble = post_preamble['n']
# debug			print "n_preamble:   "+str(n_preamble)
#		for ind_rhs,preamble in rhs_preamble.iteritems():
		for ind_rhs in range(1,n_preamble+1): 
			if tipo == 'rhs_k':
				preamble = rhs_preamble[str(ind_rhs)]
			if tipo == 'ini':
				preamble = ini_preamble[str(ind_rhs)]
			if tipo == 'post':
				preamble = post_preamble[str(ind_rhs)]
			ow(str("! "+str(ind_rhs)+" preparation phase. "+ preamble['type'] +"\n"))
			if preamble['type'] == 'compact':
				com_lhs = preamble['lhs']
				com_oper = preamble['equation'][2:4]
				com_sp_var = preamble['equation'][4]
#				index_com_sp_var = find_index_space_var(com_sp_var, space_variables)
				com_rhs = preamble['equation'][7:10]

				if tipo == 'post' or tipo == 'ini':
					com_lhs = transform_variables_per_bc(com_lhs,space_variables,variables,add_variables,sca_variables,'')
					com_rhs = transform_variables_per_bc(com_rhs,space_variables,variables,add_variables,sca_variables,'')
				else:
					com_lhs = transform_variables(com_lhs,space_variables,variables,add_variables,sca_variables)
					com_rhs = transform_variables(com_rhs,space_variables,variables,add_variables,sca_variables)

#				n_equ = variables['n']
#				for ind_var in range(1,n_equ+1):
#					name_var = variables[str(ind_var)]['symbol']
#					com_rhs = com_rhs.replace(name_var,"var_"+name_var+"_k"+matrix)
#					com_lhs = com_lhs.replace(name_var,"var_"+name_var+"_k"+matrix)
#				n_add_var = add_variables['n']
#				for ind_add_var in range(1,n_add_var+1):
#					name_add_var = add_variables[str(ind_add_var)]['symbol']
#					com_rhs = com_rhs.replace(name_add_var,"add_"+name_add_var+"_k"+matrix)
#					com_lhs = com_lhs.replace(name_add_var,"add_"+name_add_var+"_k"+matrix)
# debug				print "ehi -----------> "+com_lhs+" ; "+com_oper+" ; "+com_sp_var+" ; "+com_rhs+"\n"

				for indice in range(n_dim,0,-1):
					symbol = space_variables[str(indice)]['symbol']
					if symbol != com_sp_var:
						ow("do i_"+symbol+"=istart_"+symbol+"-ng_"+symbol+",iend_"+symbol+"+ng_"+symbol+"\n")
				stampa_compact_inversion(output,compact_operators,com_oper,com_lhs,com_rhs,space_variables,variables,add_variables,sca_variables,com_sp_var)
				for i in range(1,n_dim):  # un enddo di mesno
					ow("enddo\n")
			elif  preamble['type'] == 'internal' or preamble['type'] == 'ghost':
				ow("!PARALLEL_STANDARD\n")
				for indice in range(n_dim,0,-1):
					symbol = space_variables[str(indice)]['symbol']
					if preamble['type'] == 'internal': #tipo == 'rhs_k':
#RIMETTERE FORSE						ow("do i_"+symbol+"=istart_"+symbol+",iend_"+symbol+"\n")
						ow("do i_"+symbol+"=istart_"+symbol+",iend_"+symbol+"\n")
					if preamble['type'] == 'ghost': #tipo == 'ini' or tipo == 'post':
						ow("do i_"+symbol+"=istart_"+symbol+"-ng_"+symbol+",iend_"+symbol+"+ng_"+symbol+"\n")
# debug				print "preamble: \n "
# debug				print preamble
				n_equ_preamble = preamble['n']
# debug				print "ini_preamble: ind_rhs---> "+ str(ind_rhs)+" n: "+str(n_equ_preamble)
 				for ind_preamble in range(1,n_equ_preamble+1):
					equ_preamble = preamble[str(ind_preamble)]
					lhs = equ_preamble['lhs']
					rhs = equ_preamble['equation']
					if tipo == 'ini' or tipo == 'post':
						k_or = ''
					else:
						k_or = '_k'
					lhs = transform_variables_per_bc(lhs,space_variables,variables,add_variables,sca_variables,k_or)
					rhs = transform_variables_per_bc(rhs,space_variables,variables,add_variables,sca_variables,k_or)
					rhs = transform_expression(rhs,operators,variables,add_variables,sca_variables,space_variables)
 					#ow(indent+lhs+" = "+rhs_make(operators,variables,add_variables,sca_variables,space_variables,rhs)+"\n")
 					# random process
					if re.search("RANDOM",rhs):
						symb_1 = space_variables[str(n_dim)]['symbol']
						if n_dim > 1:
							symb_2 = space_variables[str(n_dim-1)]['symbol']
						else:
							symb_2 = ''
						stampa_random(output,symb_1,symb_2)
						rhs = re.sub('RANDOM','ran_ran2',rhs)
 					ow(indent+lhs+" = "+rhs+"\n")
#					
				for indice in range(1,n_dim+1):
					ow("enddo\n")
			elif  preamble['type'] == 'update':
# manteniamo l'update di una variabile alla volta in vista di MPI
				if tipo == 'ini' or tipo == 'post':
					k_or = ''
				else:
					k_or = '_k'
				var_update = preamble['var']
				stampa_update(output,space_variables,variables,add_variables,sca_variables,boundary,boundary_kinds,operators,k_or,var_update)
	
			elif preamble['type'] == 'scalar':
				n_equ_preamble = preamble['n']
				ow("!PARALLEL_ONLYONE_START\n")
 				for ind_preamble in range(1,n_equ_preamble+1):
					equ_preamble = preamble[str(ind_preamble)]
					lhs = equ_preamble['lhs']
					rhs = equ_preamble['equation']
					lhs = transform_variables(lhs,space_variables,variables,add_variables,sca_variables)
					rhs = transform_variables(rhs,space_variables,variables,add_variables,sca_variables)
 					# random process
					if re.search("RANDOM",rhs):
						stampa_random(output,'','')
						rhs = re.sub('RANDOM','ran_ran2',rhs)
 					ow(indent+lhs+" = "+rhs+"\n")
				ow("!PARALLEL_ONLYONE_END\n")
			elif preamble['type'] == 'sum_scalar_reduction':
				var_reduction = preamble['var']
				ow("call mpi_allreduce(MPI_IN_PLACE,sca_"+var_reduction+"_k,1,mpi_prec,MPI_SUM,MPI_COMM_WORLD, ierr)\n")
			elif preamble['type'] == 'max_scalar_reduction':
				var_reduction = preamble['var']
				ow("call mpi_allreduce(MPI_IN_PLACE,sca_"+var_reduction+"_k,1,mpi_prec,MPI_MAX,MPI_COMM_WORLD, ierr)\n")
			elif preamble['type'] == 'min_scalar_reduction':
				var_reduction = preamble['var']
				ow("call mpi_allreduce(MPI_IN_PLACE,sca_"+var_reduction+"_k,1,mpi_prec,MPI_MIN,MPI_COMM_WORLD, ierr)\n")

	if tipo != 'ini' and tipo != 'post':
		ow(str("! evaluation cycle. \n"))
		ow("!PARALLEL_STANDARD\n")
		for indice in range(n_dim,0,-1):
			symbol = space_variables[str(indice)]['symbol']
			ow("do i_"+symbol+"=istart_"+symbol+",iend_"+symbol+"\n")
		n_equ = variables['n']
		for indice in range(1,n_equ+1):
			var_symb = variables[str(indice)]['symbol']
			equation = equations[str(indice)]['equation']
	                if tipo == 'var':
				ow(indent+"var_"+var_symb+matrix+" = var_"+var_symb+matrix+" + rhs_"+var_symb+matrix+"\n")
	                if tipo == 'var_k':
				ow(indent+"var_"+var_symb+"_k"+matrix+" = var_"+var_symb+matrix+" + rhs_"+var_symb+"_k"+matrix+"*alpha(ik)*dt"+"\n")
	                if tipo == 'rhs':
				ow(indent+"rhs_"+var_symb+matrix+" = rhs_"+var_symb+matrix+" + rhs_"+var_symb+"_k"+matrix+"*beta(ik)*dt"+"\n")
	                if tipo == 'rhs_k':
				rhs = rhs_make(operators,variables,add_variables,sca_variables,space_variables,equation)
 				# random process
				if re.search("RANDOM",rhs):
					symb_1 = space_variables[str(n_dim)]['symbol']
					if n_dim > 1:
						symb_2 = space_variables[str(n_dim-1)]['symbol']
					else:
						symb_2 = ''
					stampa_random(output,symb_1,symb_2)
					rhs = re.sub('RANDOM','ran_ran2',rhs)
				ow(indent+"rhs_"+var_symb+"_k"+matrix+" = "+rhs+"\n")
		for indice in range(1,n_dim+1):
			ow("enddo\n")
	output.write("end subroutine "+tipo+"_eval\n")

def stampa_bc(output,tipo,space_variables,variables,add_variables,sca_variables,boundary,boundary_kinds,operators,mpi_or_not,block_number=''):
	k_or = tipo[3:5]
	ow("subroutine bc_"+tipo+"()\n")
	ow("use storage\n")
	if mpi_or_not == 'mpi_yes':
		ow("use mpi_storage\n")
	ow('implicit none \n')
	indici = crea_indici(space_variables)
	ow("integer "+str(indici)+"\n")
	n_dim = space_variables['n']
	n_equ = variables['n']
        for i_dim in range(1,n_dim+1):
		for start_or_end in ('start','end'):
			bou = boundary[str(i_dim)+'-'+start_or_end]
			rhs_or_var = tipo[0:3]
			bou_var = bou[rhs_or_var]
			n_bou = bou_var['n']
# debug			print "rhs_or_var,n_bou:  "+rhs_or_var+"  "+str(n_bou)+"\n"
			if n_bou != 0:
				for i_bou in range(1,n_bou+1):
					bou_i = bou_var[str(i_bou)]
					bou_i_type = bou_i['type']
					if bou_i_type == 'manual':
						bou_i_node = bou_i['node']	
						if mpi_or_not == 'mpi_yes':
							ow("if( n_truebc_"+start_or_end+"("+str(i_dim)+") .eq. 1) then\n")

						ow("!PARALLEL_REORDER\n")
						stampa_do_without_i_dim(output,space_variables,i_dim,start_or_end,bou_i_node)	
						bou_i_n = bou_i['n']
						for i_bou_i in range(1,bou_i_n+1):
							lhs = bou_i[str(i_bou_i)]['lhs']
							rhs = bou_i[str(i_bou_i)]['equation']
							lhs = transform_variables_per_bc(lhs,space_variables,variables,add_variables,sca_variables,k_or)
							rhs = transform_variables_per_bc(rhs,space_variables,variables,add_variables,sca_variables,k_or,'exchange_yes',block_number,i_dim)
							rhs = transform_expression(rhs,operators,variables,add_variables,sca_variables,space_variables)
 							if rhs_or_var == tipo[0:3]:
 								# random process
								if re.search("RANDOM",rhs):
									(symb_1,symb_2) = select_bc_random(i_dim,n_dim,space_variables)
									stampa_random(output,symb_1,symb_2)
									rhs = re.sub('RANDOM','ran_ran2',rhs)
								ow(lhs+" = "+rhs+"\n")
						stampa_enddo(output,space_variables)
						if mpi_or_not == 'mpi_yes':
							ow("endif\n")
	ow("end subroutine bc_"+tipo)

# stampa_update is just a kind of simplified stampa_bc
def stampa_update(output,space_variables,variables,add_variables,sca_variables,boundary,boundary_kinds,operators,k_or,update_var):
	n_dim = space_variables['n']
	n_equ = variables['n']
        for i_dim in range(1,n_dim+1):
		symbol = space_variables[str(i_dim)]['symbol']
		for start_or_end in ('start','end'):
			rhs_or_var = 'var'
			bou_i_node = 'ghost'
			ow("if( n_truebc_"+start_or_end+"("+str(i_dim)+") .eq. 1) then\n")
			ow("!PARALLEL_REORDER\n")
			stampa_do_without_i_dim(output,space_variables,i_dim,start_or_end,bou_i_node)	
		#update all	for i_bou_i in range(1,n_equ+1):
		#update all	lhs = variables[str(i_bou_i)]['symbol']
			lhs = update_var
			if start_or_end == 'start':
				rhs = '#[ec'+symbol+']{'+lhs+'}'
			if start_or_end == 'end':
				rhs = '#[ed'+symbol+']{'+lhs+'}'
			lhs = transform_variables_per_bc(lhs,space_variables,variables,add_variables,sca_variables,k_or)
			rhs = transform_variables_per_bc(rhs,space_variables,variables,add_variables,sca_variables,k_or)
			rhs = transform_expression(rhs,operators,variables,add_variables,sca_variables,space_variables)
			ow(lhs+" = "+rhs+"\n")
			stampa_enddo(output,space_variables)
			ow("endif\n")
	stampa_mpi_update_new(output,space_variables,variables,add_variables,sca_variables,my_kind,update_var,n_decomp)

def stampa_print_screen(output,space_variables,variables,solver_or_post):
	ow("subroutine print_screen \n")
	ow("use storage \n")
	ow("use mpi_storage \n")
	ow("implicit none \n")
	indici = crea_indici(space_variables)
	indici_half = "("+crea_indici_half(space_variables)+")"
	indici_half_per_update = "("+crea_indici_half_per_update(space_variables)+")"
	ow("integer "+str(indici)+"\n")
	n_sp_var = space_variables['n']
	first_var = variables['1']['symbol']
	ow("if( &\n")
	for ind_sp_var in range(1,n_sp_var+1):
		name_sp_var = space_variables[str(ind_sp_var)]['symbol']
		ow("(nt_"+name_sp_var+"/2.ge.istart_"+name_sp_var+").and.(nt_"+name_sp_var+"/2.le.iend_"+name_sp_var+").and. &\n")
	ow(".true.) then \n")
        ow("if(GPU) then\n")
        ow("!$acc update host(var_"+first_var+indici_half_per_update+")\n")
        ow("endif\n")
	if solver_or_post == 'solver':
		ow("     write(123,*) it,var_"+first_var+indici_half+"\n")
	if solver_or_post == 'post':
		ow("     write(123,*) post_it,var_"+first_var+indici_half+"\n")
	ow("     flush(123)"+"\n")
	ow("     write(*,*) 'var_"+first_var+indici_half+" = ', var_"+first_var+indici_half+"\n")
	ow("endif\n")
	ow("end subroutine print_screen \n")

def stampa_update_gpu(output,space_variables,variables,add_variables):
	ow("subroutine update_gpu \n")
	ow("use storage \n")
	ow("use mpi_storage \n")
	n_sp_var = space_variables['n']
	for ind_sp_var in range(1,n_sp_var+1):
		name_sp_var = space_variables[str(ind_sp_var)]['symbol']
		ow("!$acc update device(gri_g_"+name_sp_var+")\n")
	
	n_equ = variables['n']
	for ind_var in range(1,n_equ+1):
		name_var = variables[str(ind_var)]['symbol']
		ow("!$acc update device(var_"+name_var+")\n")
		ow("!$acc update device(rhs_"+name_var+")\n")
		ow("!$acc update device(var_"+name_var+"_k)\n")
		ow("!$acc update device(rhs_"+name_var+"_k)\n")
	
	n_add_var = add_variables['n']
	for ind_add_var in range(1,n_add_var+1):
		name_add_var = add_variables[str(ind_add_var)]['symbol']
		ow("!$acc update device(add_"+name_add_var+"_k)\n")
	ow("end subroutine update_gpu\n")

def stampa_update_cpu(output,space_variables,variables,add_variables):
	ow("subroutine update_cpu \n")
	ow("use storage \n")
	ow("use mpi_storage \n")

	n_equ = variables['n']
	for ind_var in range(1,n_equ+1):
		name_var = variables[str(ind_var)]['symbol']
		ow("!$acc update host(var_"+name_var+")\n")
	ow("end subroutine update_cpu\n")

def stampa_mpi_update_new(output,space_variables,variables,add_variables,sca_variables,my_kind,k_or_not_k_or_var,n_decomp):
        if n_decomp > 3:
                print "Error! MPI not available if N_decomp (=N_dim for now) greater than 3 \n"
                exit(1)
	n_equ = variables['n']
	if k_or_not_k_or_var == '' or k_or_not_k_or_var == '_k' :
		n_send = str(n_equ)
		ow("subroutine update_var"+k_or_not_k_or_var+"_mpi()\n")
		ow("\n")
		ow("use mpi\n")
		ow("use storage\n")
		ow("use mpi_storage\n")
		ow("\n")
		ow("integer :: i_x\n")
		if n_decomp >= 2:
			ow("integer :: i_y\n")
		if n_decomp == 3:
			ow("integer :: i_z\n")
	else:
		n_send = "1"
	if n_decomp == 1:
		ow("ind_buf_le = "+n_send+"*ng_x\n")
		ow("ind_buf_ri = "+n_send+"*ng_x\n")
	if n_decomp == 2:
		ow("ind_buf_le = "+n_send+"*(iend_y-istart_y+1)*ng_x\n")
		ow("ind_buf_ri = "+n_send+"*(iend_y-istart_y+1)*ng_x\n")
		ow("ind_buf_to = "+n_send+"*(iend_x-istart_x+1)*ng_y\n")
		ow("ind_buf_bo = "+n_send+"*(iend_x-istart_x+1)*ng_y\n")
	if n_decomp == 3:
		ow("ind_buf_le = "+n_send+"*(iend_y-istart_y+1)*(iend_z-istart_z+1)*ng_x\n")
		ow("ind_buf_ri = "+n_send+"*(iend_y-istart_y+1)*(iend_z-istart_z+1)*ng_x\n")
		ow("ind_buf_to = "+n_send+"*(iend_x-istart_x+1)*(iend_z-istart_z+1)*ng_y\n")
		ow("ind_buf_bo = "+n_send+"*(iend_x-istart_x+1)*(iend_z-istart_z+1)*ng_y\n")
		ow("ind_buf_fr = "+n_send+"*(iend_x-istart_x+1)*(iend_y-istart_y+1)*ng_z\n")
		ow("ind_buf_ba = "+n_send+"*(iend_x-istart_x+1)*(iend_y-istart_y+1)*ng_z\n")
	if n_decomp == 1:
		for direc in ['left','right']:
			if direc == 'left':
				ow("do i_x=istart_x,istart_x+ng_x-1\n")
				indexes = 'i_x-istart_x+1,'
			if direc == 'right':
				ow("do i_x=iend_x-ng_x+1,iend_x\n")
				indexes = 'i_x-iend_x+ng_x,'
			if k_or_not_k_or_var == '' or k_or_not_k_or_var == '_k' :
				for ind_var in range(1,n_equ+1):
					name_var = variables[str(ind_var)]['symbol']
					ow(direc+"_s_buf("+indexes+str(ind_var)+")    = var_"+name_var+k_or_not_k_or_var+"(i_x)\n")
			else:
				ow(direc+"_s_buf("+indexes+"1)    = add_"+k_or_not_k_or_var+"_k(i_x)\n")
			ow("enddo\n")
	if n_decomp == 2:
		for direc in ['left','right','bottom','top']:
			ow("!PARALLEL_REORDER\n")
			if direc == 'left':
				ow("do i_y=istart_y,iend_y\n")
				ow("do i_x=istart_x,istart_x+ng_x-1\n")
				indexes = 'i_x-istart_x+1,i_y-istart_y+1,'
			if direc == 'right':
				ow("do i_y=istart_y,iend_y\n")
				ow("do i_x=iend_x-ng_x+1,iend_x\n")
				indexes = 'i_x-iend_x+ng_x,i_y-istart_y+1,'
			if direc == 'bottom':
				ow("do i_x=istart_x,iend_x\n")
				ow("do i_y=istart_y,istart_y+ng_y-1\n")
				indexes = 'i_x-istart_x+1,i_y-istart_y+1,'
			if direc == 'top':
				ow("do i_x=istart_x,iend_x\n")
				ow("do i_y=iend_y-ng_y+1,iend_y\n")
				indexes = 'i_x-istart_x+1,i_y-iend_y+ng_y,'
			if k_or_not_k_or_var == '' or k_or_not_k_or_var == '_k' :
				for ind_var in range(1,n_equ+1):
					name_var = variables[str(ind_var)]['symbol']
					ow(direc+"_s_buf("+indexes+str(ind_var)+")    = var_"+name_var+k_or_not_k_or_var+"(i_x,i_y)\n")
			else:
				ow(direc+"_s_buf("+indexes+"1)    = add_"+k_or_not_k_or_var+"_k(i_x,i_y)\n")
			ow("enddo\n")
			ow("enddo\n")
	if n_decomp == 3:
		for direc in ['left','right','bottom','top','back','front']:
			ow("!PARALLEL_REORDER\n")
			if direc == 'left':
				ow("do i_z=istart_z,iend_z\n")
				ow("do i_y=istart_y,iend_y\n")
				ow("do i_x=istart_x,istart_x+ng_x-1\n")
				indexes = 'i_x-istart_x+1,i_y-istart_y+1,i_z-istart_z+1,'
			if direc == 'right':
				ow("do i_z=istart_z,iend_z\n")
				ow("do i_y=istart_y,iend_y\n")
				ow("do i_x=iend_x-ng_x+1,iend_x\n")
				indexes = 'i_x-iend_x+ng_x,i_y-istart_y+1,i_z-istart_z+1,'
			if direc == 'bottom':
				ow("do i_z=istart_z,iend_z\n")
				ow("do i_x=istart_x,iend_x\n")
				ow("do i_y=istart_y,istart_y+ng_y-1\n")
				indexes = 'i_x-istart_x+1,i_y-istart_y+1,i_z-istart_z+1,'
			if direc == 'top':
				ow("do i_z=istart_z,iend_z\n")
				ow("do i_x=istart_x,iend_x\n")
				ow("do i_y=iend_y-ng_y+1,iend_y\n")
				indexes = 'i_x-istart_x+1,i_y-iend_y+ng_y,i_z-istart_z+1,'
			if direc == 'back':
				ow("do i_y=istart_y,iend_y\n")
				ow("do i_x=istart_x,iend_x\n")
				ow("do i_z=istart_z,istart_z+ng_z-1\n")
				indexes = 'i_x-istart_x+1,i_y-istart_y+1,i_z-istart_z+1,'
			if direc == 'front':
				ow("do i_y=istart_y,iend_y\n")
				ow("do i_x=istart_x,iend_x\n")
				ow("do i_z=iend_z-ng_z+1,iend_z\n")
				indexes = 'i_x-istart_x+1,i_y-istart_y+1,i_z-iend_z+ng_z,'
			if k_or_not_k_or_var == '' or k_or_not_k_or_var == '_k' :
				for ind_var in range(1,n_equ+1):
					name_var = variables[str(ind_var)]['symbol']
					ow(direc+"_s_buf("+indexes+str(ind_var)+")    = var_"+name_var+k_or_not_k_or_var+"(i_x,i_y,i_z)\n")
			else:
				ow(direc+"_s_buf("+indexes+"1)    = add_"+k_or_not_k_or_var+"_k(i_x,i_y,i_z)\n")
			ow("enddo\n")
			ow("enddo\n")
			ow("enddo\n")
	ow("\n")
	ow("call mpi_sendrecv(left_s_buf ,ind_buf_le,mpi_prec,n_left  , 100, &\n")
	ow("                  right_r_buf,ind_buf_ri,mpi_prec,n_right,  100, &\n")
	ow("                  mpi_comm_world,istatus,ierr)\n")
	ow("call mpi_sendrecv(right_s_buf,ind_buf_ri,mpi_prec,n_right,  100, &\n")
	ow("                  left_r_buf ,ind_buf_le,mpi_prec,n_left  , 100, &\n")
	ow("                  mpi_comm_world,istatus,ierr)\n")
	if n_decomp >= 2:
		ow("call mpi_sendrecv(top_s_buf ,ind_buf_to,mpi_prec,n_top ,    100, &\n")
		ow("                  bottom_r_buf,ind_buf_bo,mpi_prec,n_bottom,100, &\n")
		ow("                  mpi_comm_world,istatus,ierr)\n")
		ow("call mpi_sendrecv(bottom_s_buf,ind_buf_bo,mpi_prec,n_bottom,100, &\n")
		ow("                  top_r_buf ,ind_buf_to,mpi_prec,n_top ,    100, &\n")
		ow("                  mpi_comm_world,istatus,ierr)\n")
	if n_decomp == 3:
		ow("call mpi_sendrecv(front_s_buf ,ind_buf_fr,mpi_prec,n_front ,100, &\n")
		ow("                  back_r_buf,ind_buf_ba,mpi_prec,n_back,    100, &\n")
		ow("                  mpi_comm_world,istatus,ierr)\n")
		ow("call mpi_sendrecv(back_s_buf,ind_buf_ba,mpi_prec,n_back,    100, &\n")
		ow("                  front_r_buf ,ind_buf_fr,mpi_prec,n_front ,100, &\n")
		ow("                  mpi_comm_world,istatus,ierr)\n")

	n_equ = variables['n']
	if n_decomp == 1:
		for direc in ['left','right']:
			ow("if(n_"+direc+".ge.0) then\n")
			if direc == 'left':
				ow("do i_x=istart_x-ng_x,istart_x-1\n")
				indexes = 'i_x-istart_x+ng_x+1,'
			if direc == 'right':
				ow("do i_x=iend_x+1,iend_x+ng_x\n")
				indexes = 'i_x-iend_x,'
			if k_or_not_k_or_var == '' or k_or_not_k_or_var == '_k' :
				for ind_var in range(1,n_equ+1):
					name_var = variables[str(ind_var)]['symbol']
					ow("var_"+name_var+k_or_not_k_or_var+"(i_x) = "+direc+"_r_buf("+indexes+str(ind_var)+")\n")
			else:
				ow("add_"+k_or_not_k_or_var+"_k(i_x) = "+direc+"_r_buf("+indexes+"1)\n")

			ow("enddo\n")
			ow("endif\n")
	if n_decomp == 2:
		for direc in ['left','right','bottom','top']:
			ow("if(n_"+direc+".ge.0) then\n")
			ow("!PARALLEL_REORDER\n")
			if direc == 'left':
				ow("do i_y=istart_y,iend_y\n")
				ow("do i_x=istart_x-ng_x,istart_x-1\n")
				indexes = 'i_x-istart_x+ng_x+1,i_y-istart_y+1,'
			if direc == 'right':
				ow("do i_y=istart_y,iend_y\n")
				ow("do i_x=iend_x+1,iend_x+ng_x\n")
				indexes = 'i_x-iend_x,i_y-istart_y+1,'
			if direc == 'bottom':
				ow("do i_x=istart_x,iend_x\n")
				ow("do i_y=istart_y-ng_y,istart_y-1\n")
				indexes = 'i_x-istart_x+1,i_y-istart_y+ng_y+1,'
			if direc == 'top':
				ow("do i_x=istart_x,iend_x\n")
				ow("do i_y=iend_y+1,iend_y+ng_y\n")
				indexes = 'i_x-istart_x+1,i_y-iend_y,'
			if k_or_not_k_or_var == '' or k_or_not_k_or_var == '_k' :
				for ind_var in range(1,n_equ+1):
					name_var = variables[str(ind_var)]['symbol']
					ow("var_"+name_var+k_or_not_k_or_var+"(i_x,i_y) = "+direc+"_r_buf("+indexes+str(ind_var)+")\n")
			else:
				ow("add_"+k_or_not_k_or_var+"_k(i_x,i_y) = "+direc+"_r_buf("+indexes+"1)\n")

			ow("enddo\n")
			ow("enddo\n")
			ow("endif\n")
	if n_decomp == 3:
		for direc in ['left','right','bottom','top','back','front']:
			ow("if(n_"+direc+".ge.0) then\n")
			ow("!PARALLEL_REORDER\n")
			if direc == 'left':
				ow("do i_z=istart_z,iend_z\n")
				ow("do i_y=istart_y,iend_y\n")
				ow("do i_x=istart_x-ng_x,istart_x-1\n")
				indexes = 'i_x-istart_x+ng_x+1,i_y-istart_y+1,i_z-istart_z+1,'
			if direc == 'right':
				ow("do i_z=istart_z,iend_z\n")
				ow("do i_y=istart_y,iend_y\n")
				ow("do i_x=iend_x+1,iend_x+ng_x\n")
				indexes = 'i_x-iend_x,i_y-istart_y+1,i_z-istart_z+1,'
			if direc == 'bottom':
				ow("do i_z=istart_z,iend_z\n")
				ow("do i_x=istart_x,iend_x\n")
				ow("do i_y=istart_y-ng_y,istart_y-1\n")
				indexes = 'i_x-istart_x+1,i_y-istart_y+ng_y+1,i_z-istart_z+1,'
			if direc == 'top':
				ow("do i_z=istart_z,iend_z\n")
				ow("do i_x=istart_x,iend_x\n")
				ow("do i_y=iend_y+1,iend_y+ng_y\n")
				indexes = 'i_x-istart_x+1,i_y-iend_y,i_z-istart_z+1,'
			if direc == 'back':
				ow("do i_y=istart_y,iend_y\n")
				ow("do i_x=istart_x,iend_x\n")
				ow("do i_z=istart_z-ng_z,istart_z-1\n")
				indexes = 'i_x-istart_x+1,i_y-istart_y+1,i_z-istart_z+ng_z+1,'
			if direc == 'front':
				ow("do i_y=istart_y,iend_y\n")
				ow("do i_x=istart_x,iend_x\n")
				ow("do i_z=iend_z+1,iend_z+ng_z\n")
				indexes = 'i_x-istart_x+1,i_y-istart_y+1,i_z-iend_z,'

			if k_or_not_k_or_var == '' or k_or_not_k_or_var == '_k' :
				for ind_var in range(1,n_equ+1):
					name_var = variables[str(ind_var)]['symbol']
					ow("var_"+name_var+k_or_not_k_or_var+"(i_x,i_y,i_z) = "+direc+"_r_buf("+indexes+str(ind_var)+")\n")
			else:
				ow("add_"+k_or_not_k_or_var+"_k(i_x,i_y,i_z) = "+direc+"_r_buf("+indexes+"1)\n")
			ow("enddo\n")
			ow("enddo\n")
			ow("enddo\n")
			ow("endif\n")

	if k_or_not_k_or_var == '' or k_or_not_k_or_var == '_k' :
		ow("end subroutine update_var"+k_or_not_k_or_var+"_mpi\n")

def stampa_storage_exchange(output,space_variables,variables,add_variables,block_number,exchange_data):

	ow("module storage_exchange_"+str(block_number)+"\n\n")
	ow("use basic_var \n")
	ow("use mpi_basic_var \n\n")

	n_dim = space_variables['n']
	if exchange_data['var_send'] == 'yes':
		n_exc_var_send = variables['n']
	else:
		n_exc_var_send = 0
	if exchange_data['rhs_send'] == 'yes':
		n_exc_rhs_send = variables['n']
	else:
		n_exc_rhs_send = 0
	n_exc_var_recv = exchange_data['n_var_recv']
	n_exc_rhs_recv = exchange_data['n_rhs_recv']

	index_send = 0
	if n_exc_var_send > 0 or n_exc_rhs_send > 0:
		space_var_send = exchange_data['space_var']
		index_send = find_index_space_var(space_var_send, space_variables)
	index_recv = 0
	if n_exc_var_recv > 0 or n_exc_rhs_recv > 0:
		space_var_recv = exchange_data['space_var']
		index_recv = find_index_space_var(space_var_recv, space_variables)

	mpi_proc_send = 'n_exc_send'
	mpi_var_send  = 'exc_var_s_buf_'+str(block_number)
	mpi_rhs_send  = 'exc_rhs_s_buf_'+str(block_number)
	mpi_proc_recv = 'MPI_ANY_SOURCE'
	mpi_var_recv  = 'exc_var_r_buf_'+str(block_number)
	mpi_rhs_recv  = 'exc_rhs_r_buf_'+str(block_number)

	matrix_var = "(:"
	for indice in range(1,n_dim+1):
		matrix_var = matrix_var + ",:" 
	matrix_var += ")"

	if n_exc_var_send > 0:
		ow("real(myk), dimension"+str(matrix_var)+", allocatable :: ")
		ow(mpi_var_send+"\n")
		ow("!$acc mirror("+mpi_var_send+")\n")
	if n_exc_rhs_send > 0:
		ow("real(myk), dimension"+str(matrix_var)+", allocatable :: ")
		ow(mpi_rhs_send+"\n")
		ow("!$acc mirror("+mpi_rhs_send+")\n")
	if n_exc_var_recv > 0:
		ow("real(myk), dimension"+str(matrix_var)+", allocatable :: ")
		ow(mpi_var_recv+"\n")
		ow("!$acc mirror("+mpi_var_recv+")\n")
	if n_exc_rhs_recv > 0:
		ow("real(myk), dimension"+str(matrix_var)+", allocatable :: ")
		ow(mpi_rhs_recv+"\n")
		ow("!$acc mirror("+mpi_rhs_recv+")\n")

	ow("\nend module storage_exchange_"+str(block_number)+"\n\n")

def stampa_allocate_exchange(output,space_variables,variables,add_variables,block_number,exchange_data):

	n_dim = space_variables['n']
	if exchange_data['var_send'] == 'yes':
		n_exc_var_send = variables['n']
	else:
		n_exc_var_send = 0
	if exchange_data['rhs_send'] == 'yes':
		n_exc_rhs_send = variables['n']
	else:
		n_exc_rhs_send = 0
	n_exc_var_recv = exchange_data['n_var_recv']
	n_exc_rhs_recv = exchange_data['n_rhs_recv']

	index_send = 0
	if n_exc_var_send > 0 or n_exc_rhs_send > 0:
		space_var_send = exchange_data['space_var']
		index_send = find_index_space_var(space_var_send, space_variables)
	index_recv = 0
	if n_exc_var_recv > 0 or n_exc_rhs_recv > 0:
		space_var_recv = exchange_data['space_var']
		index_recv = find_index_space_var(space_var_recv, space_variables)

	mpi_proc_send = 'n_exc_send'
	mpi_var_send  = 'exc_var_s_buf_'+str(block_number)
	mpi_rhs_send  = 'exc_rhs_s_buf_'+str(block_number)
	mpi_proc_recv = 'MPI_ANY_SOURCE'
	mpi_var_recv  = 'exc_var_r_buf_'+str(block_number)
	mpi_rhs_recv  = 'exc_rhs_r_buf_'+str(block_number)

	ow("subroutine allocate_storage_exchange_"+str(block_number)+"()\n\n")
	ow("use storage_"+str(block_number)+"\n")
	ow("use mpi_storage_"+str(block_number)+"\n")
	ow("use storage_exchange\n")
	ow("implicit none\n\n")

	indici_var_send = "("+str(n_exc_var_send)+","
	indici_var_recv = "("+str(n_exc_var_recv)+","
	indici_rhs_send = "("+str(n_exc_rhs_send)+","
	indici_rhs_recv = "("+str(n_exc_rhs_recv)+","
	for ind in range(1,n_dim+1):
		symbol = space_variables[str(ind)]['symbol']
		if ind == index_send:
			indici_var_send += 'ng_'+symbol+"+1,"
			indici_rhs_send += '1,'
		else:
			indici_var_send += 'iend_'+symbol+"-istart_"+symbol+"+1,"
			indici_rhs_send += 'iend_'+symbol+"-istart_"+symbol+"+1,"
		if ind == index_recv:
			indici_var_recv += 'ng_'+symbol+"+1,"
			indici_rhs_recv += '1,'
		else:
			indici_var_recv += 'iend_'+symbol+"-istart_"+symbol+"+1,"
			indici_rhs_recv += 'iend_'+symbol+"-istart_"+symbol+"+1,"
	indici_var_send = indici_var_send[0:len(indici_var_send)-1]+")"
	indici_var_recv = indici_var_recv[0:len(indici_var_recv)-1]+")"
	indici_rhs_send = indici_rhs_send[0:len(indici_rhs_send)-1]+")"
	indici_rhs_recv = indici_rhs_recv[0:len(indici_rhs_recv)-1]+")"

	if n_exc_var_send > 0:
		ow("allocate("+mpi_var_send+indici_var_send+")\n")
		ow(mpi_var_send+" = 0.\n")
	if n_exc_var_recv > 0:
		ow("allocate("+mpi_var_recv+indici_var_recv+")\n")
		ow(mpi_var_recv+" = 0.\n")
	if n_exc_rhs_send > 0:
		ow("allocate("+mpi_rhs_send+indici_rhs_send+")\n")
		ow(mpi_rhs_send+" = 0.\n")
	if n_exc_rhs_recv > 0:
		ow("allocate("+mpi_rhs_recv+indici_rhs_recv+")\n")
		ow(mpi_rhs_recv+" = 0.\n")

	ow("\nend subroutine allocate_storage_exchange_"+str(block_number)+"\n\n")

def stampa_mpi_exchange_get_put(output,space_variables,variables,add_variables,sca_variables,my_kind,exchange_data,n_decomp,block_number):

	n_dim = space_variables['n']
	if exchange_data['var_send'] == 'yes':
		n_exc_var_send = variables['n']
	else:
		n_exc_var_send = 0
	if exchange_data['rhs_send'] == 'yes':
		n_exc_rhs_send = variables['n']
	else:
		n_exc_rhs_send = 0
	n_exc_var_recv = exchange_data['n_var_recv']
	n_exc_rhs_recv = exchange_data['n_rhs_recv']

	index_send = 0
	if n_exc_var_send > 0 or n_exc_rhs_send > 0:
		space_var_send = exchange_data['space_var']
		index_send = find_index_space_var(space_var_send, space_variables)
	index_recv = 0
	if n_exc_var_recv > 0 or n_exc_rhs_recv > 0:
		space_var_recv = exchange_data['space_var']  # the same of send!
		index_recv = find_index_space_var(space_var_recv, space_variables)

	mpi_var_send  = 'exc_var_s_buf_'+str(block_number)
	mpi_rhs_send  = 'exc_rhs_s_buf_'+str(block_number)
	mpi_proc_recv = 'MPI_ANY_SOURCE'
	mpi_var_recv  = 'exc_var_r_buf_'+str(block_number)
	mpi_rhs_recv  = 'exc_rhs_r_buf_'+str(block_number)

       	point_index_send = exchange_data['point_index_send']

	ow_s("subroutine exc_bc_var_k_"+str(block_number)+"()\n")
	ow_s("\n")
	ow_s("use storage_"+str(block_number)+"\n")
	ow_s("use mpi_storage_"+str(block_number)+"\n")
	ow_s("use storage_exchange\n")
	ow_s("\n")
	ow_s("integer :: i_x\n")
	if n_decomp >= 2:
		ow_s("integer :: i_y\n")
	if n_decomp == 3:
		ow_s("integer :: i_z\n")
	ow_s("integer :: ind_buf_var_send, ind_buf_var_recv \n")
	ow_s("integer :: ind_buf_rhs_send, ind_buf_rhs_recv \n")
	ow_s("integer :: mpi_proc_send, point_get \n")

	size_exc_var_send = "("+str(n_exc_var_send)+" * "
	size_exc_var_recv = "("+str(n_exc_var_recv)+" * "
	size_exc_rhs_send = "("+str(n_exc_rhs_send)+" * "
	size_exc_rhs_recv = "("+str(n_exc_rhs_recv)+" * "
	mpi_proc_send = "cart_to_rank("
	for ind in range(1,n_dim+1):
		symbol = space_variables[str(ind)]['symbol']
		if ind == index_send:
			size_exc_var_send += '(ng_'+symbol+"+1)*"
			size_exc_rhs_send += '1*'
			if block_number == '001':
				mpi_proc_send += "0,"
			if block_number == '002':
				mpi_proc_send += "n_block("+str(ind)+")-1,"
		else:
			size_exc_var_send += '(iend_'+symbol+"-istart_"+symbol+"+1)*"
			size_exc_rhs_send += '(iend_'+symbol+"-istart_"+symbol+"+1)*"
			mpi_proc_send += "n_cart("+str(ind)+"),"
		if ind == index_recv:
			size_exc_var_recv += '(ng_'+symbol+"+1)*"
			size_exc_rhs_recv += '1*'
		else:
			size_exc_var_recv += '(iend_'+symbol+"-istart_"+symbol+"+1)*"
			size_exc_rhs_recv += '(iend_'+symbol+"-istart_"+symbol+"+1)*"
        size_exc_var_send = size_exc_var_send[0:len(size_exc_var_send)-1]+")"
        size_exc_var_recv = size_exc_var_recv[0:len(size_exc_var_recv)-1]+")"
        size_exc_rhs_send = size_exc_rhs_send[0:len(size_exc_rhs_send)-1]+")"
        size_exc_rhs_recv = size_exc_rhs_recv[0:len(size_exc_rhs_recv)-1]+")"
        mpi_proc_send = mpi_proc_send[0:len(mpi_proc_send)-1]+")"
	if n_exc_var_send == 0:
        	size_exc_var_send = str(0)
	if n_exc_var_recv == 0:
        	size_exc_var_recv = str(0)
	if n_exc_rhs_send == 0:
        	size_exc_rhs_send = str(0)
	if n_exc_rhs_recv == 0:
        	size_exc_rhs_recv = str(0)

	ow_s("point_get = "+str(point_index_send)+" \n")
	ow_s("mpi_proc_send = "+str(mpi_proc_send)+" \n")
	ow_s("ind_buf_var_send = "+str(size_exc_var_send)+"\n")
	ow_s("ind_buf_var_recv = "+str(size_exc_var_recv)+"\n")
	ow_s("ind_buf_rhs_send = "+str(size_exc_rhs_send)+"\n")
	ow_s("ind_buf_rhs_recv = "+str(size_exc_rhs_recv)+"\n")

	if n_exc_var_recv > 0:
		ow_s("if( n_truebc_start("+str(index_recv)+") .eq. 1) then\n")
		ow_s("call mpi_irecv("+mpi_var_recv+",ind_buf_var_recv,mpi_prec,"+mpi_proc_recv+",200, &\n")
		ow_s("                  mpi_comm_world,irequest_recv_"+str(block_number)+",ierr)\n")
		ow_s("endif\n\n")

	if n_exc_var_send > 0:
		ow_s("if((point_get.ge.istart_"+str(space_var_send)+").and.(point_get.le.iend_"+str(space_var_send)+")) then \n")
		ow_s("!PARALLEL_REORDER\n")
		stampa_do_for_exchange_mpi(output,space_variables,index_send,point_index_send,'var')
		for ind_var in range(1,n_exc_var_send+1):
	        	matrix = "("+crea_indici(space_variables)+")"
	        	matrix_buf = "("+str(ind_var)+","+ crea_indici_exchange_mpi(space_variables,'send',index_send,point_index_send)+")"
			name_var = variables[str(ind_var)]['symbol']
			ow_s("   "+mpi_var_send+matrix_buf+" = var_"+name_var+"_k"+matrix+"\n")
		stampa_enddo(output,space_variables)

		ow_s("   call mpi_isend("+mpi_var_send+",ind_buf_var_send,mpi_prec,mpi_proc_send,200, &\n")
		ow_s("                  mpi_comm_world,irequest_send_"+str(block_number)+",ierr)\n")
		ow_s("endif\n\n")

	ow_s("end subroutine exc_bc_var_k_"+str(block_number)+"\n")

	ow_s("subroutine exc_bc_rhs_k_"+str(block_number)+"()\n")
	ow_s("\n")
	ow_s("use storage_"+str(block_number)+"\n")
	ow_s("use mpi_storage_"+str(block_number)+"\n")
	ow_s("use storage_exchange\n")
	ow_s("\n")
	ow_s("integer :: i_x\n")
	if n_decomp >= 2:
		ow_s("integer :: i_y\n")
	if n_decomp == 3:
		ow_s("integer :: i_z\n")
	ow_s("integer :: ind_buf_var_send, ind_buf_var_recv \n")
	ow_s("integer :: ind_buf_rhs_send, ind_buf_rhs_recv \n")
	ow_s("integer :: mpi_proc_send, point_get \n")

	size_exc_var_send = "("+str(n_exc_var_send)+" * "
	size_exc_var_recv = "("+str(n_exc_var_recv)+" * "
	size_exc_rhs_send = "("+str(n_exc_rhs_send)+" * "
	size_exc_rhs_recv = "("+str(n_exc_rhs_recv)+" * "
	mpi_proc_send = "cart_to_rank("
	for ind in range(1,n_dim+1):
		symbol = space_variables[str(ind)]['symbol']
		if ind == index_send:
			size_exc_var_send += '(ng_'+symbol+"+1)*"
			size_exc_rhs_send += '1*'
			if block_number == '001':
				mpi_proc_send += "0,"
			if block_number == '002':
				mpi_proc_send += "n_block("+str(ind)+")-1,"
		else:
			size_exc_var_send += '(iend_'+symbol+"-istart_"+symbol+"+1)*"
			size_exc_rhs_send += '(iend_'+symbol+"-istart_"+symbol+"+1)*"
			mpi_proc_send += "n_cart("+str(ind)+"),"
		if ind == index_recv:
			size_exc_var_recv += '(ng_'+symbol+"+1)*"
			size_exc_rhs_recv += '1*'
		else:
			size_exc_var_recv += '(iend_'+symbol+"-istart_"+symbol+"+1)*"
			size_exc_rhs_recv += '(iend_'+symbol+"-istart_"+symbol+"+1)*"
        size_exc_var_send = size_exc_var_send[0:len(size_exc_var_send)-1]+")"
        size_exc_var_recv = size_exc_var_recv[0:len(size_exc_var_recv)-1]+")"
        size_exc_rhs_send = size_exc_rhs_send[0:len(size_exc_rhs_send)-1]+")"
        size_exc_rhs_recv = size_exc_rhs_recv[0:len(size_exc_rhs_recv)-1]+")"
        mpi_proc_send = mpi_proc_send[0:len(mpi_proc_send)-1]+")"
	if n_exc_var_send == 0:
        	size_exc_var_send = str(0)
	if n_exc_var_recv == 0:
        	size_exc_var_recv = str(0)
	if n_exc_rhs_send == 0:
        	size_exc_rhs_send = str(0)
	if n_exc_rhs_recv == 0:
        	size_exc_rhs_recv = str(0)

	ow_s("point_get = "+str(point_index_send)+" \n")
	ow_s("mpi_proc_send = "+str(mpi_proc_send)+" \n")
	ow_s("ind_buf_var_send = "+str(size_exc_var_send)+"\n")
	ow_s("ind_buf_var_recv = "+str(size_exc_var_recv)+"\n")
	ow_s("ind_buf_rhs_send = "+str(size_exc_rhs_send)+"\n")
	ow_s("ind_buf_rhs_recv = "+str(size_exc_rhs_recv)+"\n")

	if n_exc_rhs_recv > 0:
		ow_s("if( n_truebc_start("+str(index_recv)+") .eq. 1) then\n")
		ow_s("call mpi_irecv("+mpi_rhs_recv+",ind_buf_rhs_recv,mpi_prec,"+mpi_proc_recv+",200, &\n")
		ow_s("                  mpi_comm_world,irequest_recv_"+str(block_number)+",ierr)\n")
		ow_s("endif\n\n")

	if n_exc_rhs_send > 0:
		ow_s("if((point_get.ge.istart_"+str(space_var_send)+").and.(point_get.le.iend_"+str(space_var_send)+")) then \n")
		ow_s("!PARALLEL_REORDER\n")
		stampa_do_for_exchange_mpi(output,space_variables,index_send,point_index_send,'rhs')
		for ind_var in range(1,n_exc_rhs_send+1):
	        	matrix = "("+crea_indici(space_variables)+")"
	        	matrix_buf = "("+str(ind_var)+","+ crea_indici_exchange_mpi(space_variables,'send',index_send,point_index_send)+")"
			name_var = variables[str(ind_var)]['symbol']
			ow_s("   "+mpi_rhs_send+matrix_buf+" = rhs_"+name_var+"_k"+matrix+"\n")
		stampa_enddo(output,space_variables)

		ow_s("   call mpi_isend("+mpi_rhs_send+",ind_buf_rhs_send,mpi_prec,mpi_proc_send,200, &\n")
		ow_s("                  mpi_comm_world,irequest_send_"+str(block_number)+",ierr)\n")
		ow_s("endif\n\n")


	ow_s("end subroutine exc_bc_rhs_k_"+str(block_number)+"\n")

	ow_s("\nsubroutine exc_wait_var_k_"+str(block_number)+"\n\n")

        ow_s("use storage_"+str(block_number)+"\n")
        ow_s("use mpi_storage_"+str(block_number)+"\n")
        ow_s("use storage_exchange\n\n")

        ow_s("integer :: point_get \n\n")

	ow_s("point_get = "+str(point_index_send)+" \n")

	if n_exc_var_send > 0:
		ow_s("if((point_get.ge.istart_"+str(space_var_send)+").and.(point_get.le.iend_"+str(space_var_send)+")) then \n")
		ow_s("   call mpi_wait(irequest_send_"+str(block_number)+", istatus_send_"+str(block_number)+", ierr) ! var exchange\n")
		ow_s("endif\n\n")
	if n_exc_var_recv > 0:
		ow_s("if( n_truebc_start("+str(index_recv)+") .eq. 1) then\n")
		ow_s("   call mpi_wait(irequest_recv_"+str(block_number)+", istatus_recv_"+str(block_number)+", ierr) ! var exchange\n")
		ow_s("endif\n\n")

	ow_s("\nend subroutine exc_wait_var_k_"+str(block_number)+"\n")

	ow_s("\nsubroutine exc_wait_rhs_k_"+str(block_number)+"\n\n")

        ow_s("use storage_"+str(block_number)+"\n")
        ow_s("use mpi_storage_"+str(block_number)+"\n")
        ow_s("use storage_exchange\n\n")

        ow_s("integer :: point_get \n\n")

	ow_s("point_get = "+str(point_index_send)+" \n")

	if n_exc_rhs_send > 0:
		ow_s("if((point_get.ge.istart_"+str(space_var_send)+").and.(point_get.le.iend_"+str(space_var_send)+")) then \n")
		ow_s("   call mpi_wait(irequest_send_"+str(block_number)+", istatus_send_"+str(block_number)+", ierr) ! rhs exchange\n")
		ow_s("endif\n\n")
	if n_exc_rhs_recv > 0:
		ow_s("if( n_truebc_start("+str(index_recv)+") .eq. 1) then\n")
		ow_s("   call mpi_wait(irequest_recv_"+str(block_number)+", istatus_recv_"+str(block_number)+", ierr) ! rhs exchange\n")
		ow_s("endif\n\n")

	ow_s("\nend subroutine exc_wait_rhs_k_"+str(block_number)+"\n")

#MAIN
#--------------------------------------------------------------------------
# OLD
#--------------------------------------------------------------------------
# just for  ASCII I/O, do not take care of that
visual = 'tecplot' # 'gnuplot'
# just a list of possible predefined boundary types, not yet implemented
boundary_kinds = {
# var boundaries: ghost and/or boundary
# rhs boundaries: boundary
# mixed conditions: var[ghost]+rhs[boundary]
      'free':'none',               #   nothing
      'fixed':'rhs',               #   boundary
      'periodicity':'var',         #   ghost+boundary
      'extr_cost':'var',           #   ghost
      'extr_lin':'var',            #   ghost
      'dirichelet':'var',          #   ghost or ghost+boundary, value required
      'neumann':'var',             #   ghost or ghost+boundary, value required
      'manual-var':'var',          #   ghost or ghost+boundary
      'manual-rhs':'rhs',          #   boundary
      'robin':'var',               #   .....
      'orlansky':'rhs',            #   .....
      'nonreflecting':'rhs',       #   .....
      'tractionfree':'var',        #   .....
      'shock':'var'                #   ---> manuale
}
#--------------------------------------------------------------------------
#--------------------------------------------------------1
# TEST CASE WAS DEFINED HERE!!
#--------------------------------------------------------1
# add post_add_variables to add_variables avoiding duplications
N_add_var_more = 0
for i in range(1,post_add_variables['n']+1):
	add_i_var = 'yes'
	for k in range(1,add_variables['n']+1):
		if post_add_variables[str(i)]['symbol'] == add_variables[str(k)]['symbol']:
			add_i_var = 'no'
	if add_i_var == 'yes':
		N_add_var_more += 1
		add_variables[str(add_variables['n']+N_add_var_more)] = post_add_variables[str(i)]
add_variables['n'] = add_variables['n'] + N_add_var_more

#------------------------------------------------------------------------
#WEB_TEST_DIM_BODY
#------------------------------------------------------------------------

# -----------------------------------------------------------------------
# BASE 
# -----------------------------------------------------------------------
n_decomp = N_dim # MPI decomposition over all the spatial dimensions
all_private_list=[]

mem_used = memory_usage_calc(my_kind,N_dim,N_equ,N_add_var,space_variables)
print "The selected run requires total amount of memory: "+str(round(mem_used,3))+" Gb\n"

assegna_nodi_seriale(space_variables)

stampa_grid('fortran/grid',space_variables)

# print 'Number of arguments:', len(sys.argv), 'arguments.'
# print 'Argument List:', str(sys.argv)
n_command_line_arguments = len(sys.argv)
argument_list = sys.argv
if n_command_line_arguments == 2:
	if argument_list[1] == 'grid':
		print "Mode only grid generation. Stopping now"
		exit(1)

#exit(1)
#raw_input("Press ENTER to continue")

output = open("fortran/basic_var_acc.F90", 'w')
stampa_basic_var(output, my_kind, space_variables, 'post_no',n_decomp)
stampa_time_init(output,i_kutta)


output = open("fortran/storage_allocatable.F90", 'w')
stampa_storage_allocatable(output,i_kutta,my_kind,space_variables,variables,add_variables,sca_variables,compact_operators,n_decomp,mpi_periodic,all_private_list)
output.close()
output = open("fortran/allocate_storage.F90", 'w')
stampa_allocate_storage(output,i_kutta,my_kind,space_variables,variables,add_variables,sca_variables,compact_operators,n_decomp,random_ornot, pos_variables, average_dir)
output.close()

output = open("fortran/rk-init.f90", 'w')
stampa_size_init(output,space_variables)
stampa_random_init(output,random_ornot)
stampa_init(output,space_variables,variables,add_variables, sca_variables,initial,'var')
stampa_init(output,space_variables,variables,add_variables, sca_variables,initial,'var_k')
stampa_init(output,space_variables,variables,add_variables, sca_variables,initial,'rhs_k')
stampa_eval(output,space_variables,variables,add_variables,'ini'  ,rhs_preamble,ini_preamble,sca_variables)
output.close()

output = open("fortran/comp_init.f90", 'w')
stampa_compact_init(output,compact_operators,space_variables,variables,add_variables,sca_variables)
output.close()

output = open("fortran/rk-eval.f90", 'w')
stampa_init(output,space_variables,variables,add_variables, sca_variables,initial,'rhs')
stampa_eval(output,space_variables,variables,add_variables,'post' ,rhs_preamble,ini_preamble,sca_variables)
stampa_eval(output,space_variables,variables,add_variables,'var'  ,rhs_preamble,ini_preamble,sca_variables)
stampa_eval(output,space_variables,variables,add_variables,'var_k',rhs_preamble,ini_preamble,sca_variables)
stampa_eval(output,space_variables,variables,add_variables,'rhs'  ,rhs_preamble,ini_preamble,sca_variables)
stampa_eval(output,space_variables,variables,add_variables,'rhs_k',rhs_preamble,ini_preamble,sca_variables)
output.close()

# -----------------------------------------------------------------------
# MPI
# -----------------------------------------------------------------------
mpi = 1
# in piu' attezione a mpi_yes messo sopra
if mpi == 1:
	output = open("fortran/inpout_MPI.F90", 'w')
	stampa_read_grid(output,space_variables,'non_tot')
	stampa_read_grid(output,space_variables,'tot')
	stampa_print_sol(output,space_variables,variables,visual,'mpi_yes')
	stampa_read_sol(output,space_variables,variables,visual,'mpi_yes')
	stampa_print_sol_binary(output,space_variables,variables,sca_variables,visual,'mpi_yes')
	stampa_read_sol_binary(output,space_variables,variables,sca_variables,visual,'mpi_yes')
	stampa_writeread_sol_binary_pio(output,space_variables,variables,sca_variables,visual,'mpi_yes','write')
	stampa_writeread_sol_binary_pio(output,space_variables,variables,sca_variables,visual,'mpi_yes','read')
	output.close()
	output = open("fortran/bc_var_MPI.f90", 'w')
	stampa_bc(output,'var',space_variables,variables,add_variables,sca_variables,boundary,boundary_kinds,operators,'mpi_yes',block_number)
	output.close()
	output = open("fortran/bc_var_k_MPI.f90", 'w')
	stampa_bc(output,'var_k',space_variables,variables,add_variables,sca_variables,boundary,boundary_kinds,operators,'mpi_yes',block_number)
	output.close()
	output = open("fortran/bc_rhs_MPI.f90", 'w')
	stampa_bc(output,'rhs',space_variables,variables,add_variables,sca_variables,boundary,boundary_kinds,operators,'mpi_yes',block_number)
	output.close()
	output = open("fortran/bc_rhs_k_MPI.f90", 'w')
	stampa_bc(output,'rhs_k',space_variables,variables,add_variables,sca_variables,boundary,boundary_kinds,operators,'mpi_yes',block_number)
	output.close()

	output = open("fortran/print_screen_MPI.f90", 'w')
	stampa_print_screen(output,space_variables,variables,'solver')
	output.close()
        output = open("fortran/output_visual.F90", 'w')
	stampa_print_vtk(output, space_variables, variables, average_dir,my_kind,'solver')

# -----------------------------------------------------------------------
# GPU ACCELERATOR
# -----------------------------------------------------------------------
gpu_acc = 1

if gpu_acc == 1:
	output = open("fortran/storage_allocatable_acc.F90", 'w')
	stampa_modify_storage_per_gpu_acc_openmp(output,"fortran/storage_allocatable_acc.F90",i_kutta,space_variables,variables,add_variables,sca_variables,compact_operators,n_decomp,random_ornot,'gpuacc',all_private_list)

	output = open("fortran/update_acc.f90", 'w')
	stampa_update_gpu(output,space_variables,variables,add_variables)
	stampa_update_cpu(output,space_variables,variables,add_variables)

	output = open("fortran/bc_var_MPI_acc.f90", 'w')
	stampa_modify_per_gpu_acc_openmp(output,"fortran/bc_var_MPI_acc.f90",i_kutta,space_variables,variables,add_variables,sca_variables,compact_operators,n_decomp,'gpuacc',all_private_list)
	output = open("fortran/bc_var_k_MPI_acc.f90", 'w')
	stampa_modify_per_gpu_acc_openmp(output,"fortran/bc_var_k_MPI_acc.f90",i_kutta,space_variables,variables,add_variables,sca_variables,compact_operators,n_decomp,'gpuacc',all_private_list)
	output = open("fortran/bc_rhs_MPI_acc.f90", 'w')
	stampa_modify_per_gpu_acc_openmp(output,"fortran/bc_rhs_MPI_acc.f90",i_kutta,space_variables,variables,add_variables,sca_variables,compact_operators,n_decomp,'gpuacc',all_private_list)
	output = open("fortran/bc_rhs_k_MPI_acc.f90", 'w')
	stampa_modify_per_gpu_acc_openmp(output,"fortran/bc_rhs_k_MPI_acc.f90",i_kutta,space_variables,variables,add_variables,sca_variables,compact_operators,n_decomp,'gpuacc',all_private_list)
	output = open("fortran/bc_var_MPI_omp.f90", 'w')
	stampa_modify_per_gpu_acc_openmp(output,"fortran/bc_var_MPI_omp.f90",i_kutta,space_variables,variables,add_variables,sca_variables,compact_operators,n_decomp,'openmp',all_private_list)
	output = open("fortran/bc_var_k_MPI_omp.f90", 'w')
	stampa_modify_per_gpu_acc_openmp(output,"fortran/bc_var_k_MPI_omp.f90",i_kutta,space_variables,variables,add_variables,sca_variables,compact_operators,n_decomp,'openmp',all_private_list)
	output = open("fortran/bc_rhs_MPI_omp.f90", 'w')
	stampa_modify_per_gpu_acc_openmp(output,"fortran/bc_rhs_MPI_omp.f90",i_kutta,space_variables,variables,add_variables,sca_variables,compact_operators,n_decomp,'openmp',all_private_list)
	output = open("fortran/bc_rhs_k_MPI_omp.f90", 'w')
	stampa_modify_per_gpu_acc_openmp(output,"fortran/bc_rhs_k_MPI_omp.f90",i_kutta,space_variables,variables,add_variables,sca_variables,compact_operators,n_decomp,'openmp',all_private_list)

	output = open("fortran/rk-eval_acc.f90", 'w')
	stampa_modify_per_gpu_acc_openmp(output,"fortran/rk-eval_acc.f90",i_kutta,space_variables,variables,add_variables,sca_variables,compact_operators,n_decomp,'gpuacc',all_private_list,'','mpi_only_one')
	output = open("fortran/rk-eval_omp.f90", 'w')
	stampa_modify_per_gpu_acc_openmp(output,"fortran/rk-eval_omp.f90",i_kutta,space_variables,variables,add_variables,sca_variables,compact_operators,n_decomp,'openmp',all_private_list)

	output = open("fortran/mpi_update_new.f90", 'w')
	stampa_mpi_update_new(output,space_variables,variables,add_variables,sca_variables,my_kind,'',n_decomp)
	stampa_mpi_update_new(output,space_variables,variables,add_variables,sca_variables,my_kind,'_k',n_decomp)
	output = open("fortran/mpi_update_new_acc.f90", 'w')
	stampa_modify_per_gpu_acc_openmp(output,"fortran/mpi_update_new_acc.f90",i_kutta,space_variables,variables,add_variables,sca_variables,compact_operators,n_decomp,'gpuacc',all_private_list)
	output = open("fortran/mpi_update_new_omp.f90", 'w')
	stampa_modify_per_gpu_acc_openmp(output,"fortran/mpi_update_new_omp.f90",i_kutta,space_variables,variables,add_variables,sca_variables,compact_operators,n_decomp,'openmp',all_private_list)
	output = open("fortran/storage_allocatable_omp.F90", 'w')
	stampa_modify_storage_per_gpu_acc_openmp(output,"fortran/storage_allocatable_omp.F90",i_kutta,space_variables,variables,add_variables,sca_variables,compact_operators,n_decomp,random_ornot,'openmp',all_private_list)

# -----------------------------------------------------------------------
# POST-PROCESSING
# -----------------------------------------------------------------------
if post_proc == "yes":
	if post_proc_kind == 'avg_xyz':
		shutil.copy("fortran/post/post-fft_xyz.f90","fortran/post/post-fft.f90")
	if post_proc_kind == 'avg_yz':
		shutil.copy("fortran/post/post-fft_yz.f90","fortran/post/post-fft.f90")
	output = open("fortran/post/storage_allocatable.F90", 'w')
	stampa_storage_allocatable_post(output,i_kutta,my_kind,space_variables,variables,add_variables,sca_variables,compact_operators,n_decomp,mpi_periodic,pos_variables,tot_variables,average_dir)

	output = open("fortran/post/allocate_storage.F90", 'w')
	stampa_allocate_storage_post(output,i_kutta,my_kind,space_variables,variables,add_variables,sca_variables,compact_operators,n_decomp,random_ornot,pos_variables,tot_variables,average_dir)

	output = open("fortran/post/basic_var_acc.F90", 'w')
	stampa_basic_var(output, my_kind, space_variables, 'post_yes',n_decomp)

        file_noblock = "fortran/inpout_MPI.F90"
        file_temp = file_noblock+".temp"
        file_block = file_noblock+".block"
	shutil.copy(file_noblock,file_temp)
	output = open(file_temp,'a')
#	stampa_read_grid(output,space_variables,'tot')
	stampa_writeread_sol_binary_pio_stream(output,space_variables,variables,sca_variables,visual,'mpi_yes','read')
	stampa_writeread_sol_binary_pio_stream(output,space_variables,variables,sca_variables,visual,'mpi_yes','write')
        output.close()
	put_block_constants_post_and_post_it(file_temp,file_block)
	shutil.move(file_block,"fortran/post/inpout_MPI.F90")
	os.remove(file_temp)

	output = open("fortran/post/post-init-eval.F90", 'w')
	stampa_init_post(output,pos_variables)
	stampa_init_tot(output, tot_variables)

	stampa_post(output,space_variables,post_mpi_first,'mpi' ,'first', average_dir,variables,add_variables,sca_variables)
	stampa_post(output,space_variables,post_mpi_second,'mpi' ,'second',average_dir,variables,add_variables,sca_variables)
	stampa_post(output,space_variables,post_full_first,'full','first', average_dir,variables,add_variables,sca_variables)
	stampa_post(output,space_variables,post_full_second,'full','second',average_dir,variables,add_variables,sca_variables)

	stampa_reduce_post(output, space_variables,reduce_mpi_first,  'mpi','first'  ,average_dir,'space_and_time','reduce')
	stampa_reduce_post(output, space_variables,reduce_mpi_second, 'mpi','second' ,average_dir,'space_and_time','reduce')
	stampa_reduce_post(output, space_variables,reduce_mpi_second, 'mpi','second' ,average_dir,'space','reduce')
	stampa_reduce_post(output, space_variables,reduce_full_first, 'full','first' ,average_dir,'space_and_time','reduce')
	stampa_reduce_post(output, space_variables,reduce_full_second,'full','second',average_dir,'space_and_time','reduce')
	stampa_reduce_post(output, space_variables,reduce_full_second,'full','second',average_dir,'space','reduce')
#       for run-time post-processing
	stampa_reduce_post(output, space_variables,reduce_mpi_first, 'mpi','first' ,average_dir,'space_and_time','enlarge')
	stampa_reduce_post(output,space_variables, reduce_mpi_second,'mpi','second',average_dir,'space_and_time','enlarge')

	stampa_print_average_tecplot(output, space_variables, pos_variables, average_dir,'space_and_time')
	stampa_print_average_tecplot(output, space_variables, pos_variables, average_dir,'space')
	stampa_read_average_tecplot(output, space_variables, pos_variables, average_dir,'space_and_time')

	output.close()

	shutil.copy("fortran/post/post-init-eval.F90","fortran/post-init-eval.F90")

	output = open("fortran/post/reconstr_MPI.F90", 'w')
	stampa_size_init(output,space_variables)
	stampa_total_print_tecplot(output, space_variables, variables, average_dir)
	stampa_total_print_plot3d(output, space_variables, variables, average_dir)
	stampa_total_print_vtk(output, space_variables, variables, average_dir)
	stampa_print_vtk(output, space_variables, variables, average_dir,my_kind,'post')

	output = open("fortran/post/print_screen_MPI.f90", 'w')
	stampa_print_screen(output,space_variables,variables,'post')
	output.close()

# -----------------------------------------------------------------------
# -----------------------------------------------------------------------
solver_sources_list = ['storage_allocatable.F90','mpi_utils_new.f90','inpout_MPI.F90','rk-init.f90','comp_init.f90','rk-eval.f90','bc_var_MPI.f90','bc_var_k_MPI.f90','bc_rhs_MPI.f90','bc_rhs_k_MPI.f90','allocate_storage.F90','storage_allocatable_acc.F90','storage_allocatable_omp.F90','bc_var_MPI_acc.f90','bc_var_k_MPI_acc.f90','bc_rhs_MPI_acc.f90','bc_rhs_k_MPI_acc.f90','rk-eval_acc.f90','mpi_update_new.f90','mpi_update_new_acc.f90','bc_var_MPI_omp.f90','bc_var_k_MPI_omp.f90','bc_rhs_MPI_omp.f90','bc_rhs_k_MPI_omp.f90','rk-eval_omp.f90','mpi_update_new_omp.f90','output_visual.F90']
# -----------------------------------------------------------------------
# -----------------------------------------------------------------------

# -----------------------------------------------------------------------
# MANUAL
# -----------------------------------------------------------------------
manual_subroutines=[]
file_manual_name = "fortran/manual.F90"
if os.path.exists(file_manual_name):
	file_manual = open(file_manual_name, 'r')
        f = file_manual.readlines()
        for line in f:
                if re.search("end subroutine",line):
			manual_subroutine_name = re.sub('end subroutine ([a-z]+)','\\1',line)
			manual_subroutines.append(manual_subroutine_name.rstrip())  # strip end new lines
	file_manual.close()

	for file_withmanual_name in solver_sources_list:
		file_withmanual_name    = "fortran/"+file_withmanual_name
		file_withoutmanual_name = file_withmanual_name+".nomanual"
		file_withmanual = open(file_withmanual_name, 'r')
		file_withoutmanual = open(file_withoutmanual_name, 'w')
		delete = 0
		delete_next = 1
		for line in file_withmanual:
			for manual_subroutine_name in manual_subroutines:
				if re.search("( )*subroutine( )*"+manual_subroutine_name,line):
					delete = 1
				if re.search("( )*end\s*subroutine( )*"+manual_subroutine_name,line):
					delete_next = 0
			if delete == 0:
				file_withoutmanual.write(line)
			if delete_next == 0:
				delete_next = 1
				delete = 0
                shutil.move(file_withoutmanual_name,file_withmanual_name)
		file_withmanual.close()
		file_withoutmanual.close()

print "Manual Subroutines (manual.F90 file) :"
print manual_subroutines

# -----------------------------------------------------------------------
# SETTING KIND
# -----------------------------------------------------------------------
for file_nokind in solver_sources_list:
	output.close()
	file_nokind = "fortran/"+file_nokind
	if os.path.exists(file_nokind):
		file_kind = file_nokind+".kind"
		put_kind_constants(file_nokind,file_kind)
		shutil.move(file_kind,file_nokind)

for file_nokind in ['reconstr_MPI.F90','storage_allocatable.F90','allocate_storage.F90','basic_var_acc.F90','post-init-eval.F90','inpout_MPI.F90']:
	output.close()
	file_nokind = "fortran/post/"+file_nokind
	if os.path.exists(file_nokind):
		file_kind = file_nokind+".kind"
		put_kind_constants(file_nokind,file_kind)
		shutil.move(file_kind,file_nokind)
# -----------------------------------------------------------------------
# SETTING BLOCK
# -----------------------------------------------------------------------
if block_number != "":
	for file_noblock in ['inpout_MPI.F90','rk-init.f90','comp_init.f90','rk-eval.f90','bc_var_MPI.f90','bc_var_k_MPI.f90','bc_rhs_MPI.f90','bc_rhs_k_MPI.f90','allocate_storage.F90','storage_allocatable_acc.F90','storage_allocatable_omp.F90','bc_var_MPI_acc.f90','bc_var_k_MPI_acc.f90','bc_rhs_MPI_acc.f90','bc_rhs_k_MPI_acc.f90','rk-eval_acc.f90','mpi_update_new.f90','mpi_update_new_acc.f90','update_acc.f90','print_screen_MPI.f90','grid_x.dat','grid_y.dat','grid_z.dat','mpi_decompose.F90','bc_var_MPI_omp.f90','bc_var_k_MPI_omp.f90','bc_rhs_MPI_omp.f90','bc_rhs_k_MPI_omp.f90','rk-eval_omp.f90','mpi_update_new_omp.f90','scalars.f90','post-init-eval.F90','output_visual.F90','manual.F90']:
		output.close()
		file_noblock_np = file_noblock
		file_noblock = "fortran/"+file_noblock
		if os.path.exists(file_noblock):
			file_block = file_noblock+".block"
			put_block_constants(file_noblock,file_block,block_number)
			shutil.move(file_block,"fortran/"+block_number+"_"+file_noblock_np)

	output = open("fortran/"+str(block_number)+"_storage_exchange_acc.f90", 'w')
	stampa_storage_exchange(output,space_variables,variables,add_variables,block_number,exchange_data)
	output = open("fortran/"+str(block_number)+"_allocate_exchange.f90", 'w')
	stampa_allocate_exchange(output,space_variables,variables,add_variables,block_number,exchange_data)
	output = open("fortran/"+str(block_number)+"_exchange_get_put.f90", 'w')
	stampa_mpi_exchange_get_put(output,space_variables,variables,add_variables,sca_variables,my_kind,exchange_data,n_decomp,block_number)
	if gpu_acc == 1:
		output = open("fortran/"+str(block_number)+"_exchange_get_put_acc.f90", 'w')
		stampa_modify_per_gpu_acc_openmp(output,"fortran/"+str(block_number)+"_exchange_get_put_acc.f90",i_kutta,space_variables,variables,add_variables,sca_variables,compact_operators,n_decomp,'gpuacc',all_private_list,block_number)
		output = open("fortran/"+str(block_number)+"_exchange_get_put_omp.f90", 'w')
		stampa_modify_per_gpu_acc_openmp(output,"fortran/"+str(block_number)+"_exchange_get_put_omp.f90",i_kutta,space_variables,variables,add_variables,sca_variables,compact_operators,n_decomp,'openmp',all_private_list,block_number)

#print "all_private_list", all_private_list

exit(1)
#
